PREFACE 

The  enclosed  report  is  the  work  of  Mr.  Tao-Ming  Chang  the  graduate 
student  working  under  this  research  project.  The  report  is  entitled 
Evaluation  of  surface-wave  waveform  modeling  for  lithosphere  velocity 
structure.  The  text  is  the  draft  of  a  Ph.  D.  Dissertation. 

This  report  is  the  culmination  of  efforts  to  implement  and  compare 
broadband  waveform  modeling  techniques  to  constrain  crustal  velocity 
models.  The  impetus  for  this  effort  lies  in  the  need  of  CTBT  monitor¬ 
ing  for  accurate  event  location  and  confident  identification  from  sparse 
data.  The  velocity  models  obtained  by  analyzing  large  seismic  events 
may  provide  the  necessary  regional  constraints  required  for  the  small 
events  of  interest  to  the  CTBT. 

This  report  implements  and  compares  techniques  for  one  well  recorded 
North  American  earthquake  rather  than  looking  at  a  lot  of  earth¬ 
quakes.  There  is  enough  variability  in  wave  propagation  in  North 
America  to  test  a  method.  The  single  event  is  understood  well  enough 
that  the  effects  of  the  source  on  the  waveform  is  constrained. 

This  effort  is  just  the  beginning  of  a  long  effort  of  algorithm  testing  and 
improvement. 
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2.1  Station  distribution  map.  The  stations  are  shown  as  trian-  13 
gles.  The  focal  mechanism  used  in  this  study  is  shown  as 

the  beach  ball.  This  map  is  generated  using  GMT  (Wessel 
and  Smith,  1991) 

2.2  For  teleseismic  records,  a  general  7-9  seconds  travel  time  14 
difference  for  pP-P  is  observable. 

2.3  The  Love  wave  radiation  patterns  for  the  preferred  focal  15 
mechanism  at  23  km  depth.  The  bars  indicate  attenuation 
corrected  spectral  amplitudes  in  cm-sec  normalized  for 
geometrical  spreading  to  1000  km. 

2.4  The  Love  wave  radiation  patterns  for  Harvard  CMT  solu-  16 
tion  at  15  km  depth. 

2.5  The  Love  wave  radiation  patterns  for  Sipkin’s  USGS  solu-  17 
tion  at  13  km  depth. 

2.6  The  Rayleigh  wave  radiation  patterns  for  the  preferred  18 
focal  mechanism  at  23  km  depth. 

2.7  The  Rayleigh  wave  radiation  patterns  for  Harvard  CMT  19 
solution  at  15  km  depth. 

2.8  The  Rayleigh  wave  radiation  patterns  for  Sipkin’s  USGS  20 
solution  at  13  km  depth. 

2.9  This  figure  shows  the  observed  P-wave  polarities  and  the  21 
reestimated  focal  mechanism  which  will  be  used  for  test¬ 
ing  several  surface-waveform  modeling  algorithms. 

2.10  This  figure  shows  the  focal  mechanism  of  Harvard  CMT  22 
solution  and  the  observed  P-wave  polarities. 

2.11  This  figm-e  shows  the  focal  mechanism  of  Sipkin’s  USGS  23 
solution  and  the  observed  P-wave  polarities. 
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3.1a  This  is  an  example  showning  how  unexpected  signal  42 
interference  affects  the  extracting  procedure  in  GSDF  the¬ 
ory.  There  are  five  traces  presented  to  show  the  different 
processing  stages.  The  top  two  traces  are  the  prefiltered 
isolation  filter  and  observed  seismogram,  respectively, 

The  third  trace  shows  the  Gaussian  filtered  cross¬ 
correlation  at  a  target  period  of  20.0  seconds.  The  five 
extracted  parameters  are  shown.  The  dashed  curves 
inside  the  envelope  are  from  the  synthetics  and  the  solid 
curve  are  from  the  observed  data.  The  fourth  trace  is  the 
windowed  cross-correlagram  shown  in  the  third  trace.  The 
bottom  trace  is  the  filtered  windowed  cross-correlagram 
and  its  five  Gaussian  wavelet  parameters  which  are  to  be 
used  in  further  processing.  Due  to  improper  rotation,  the 
Love  wave  appears  on  radial  component  before  the 
Rayleigh  wave  arrival.  The  Love  wave  wavetrain  may 
causes  two  kind  troubles  like  (a)  wrong  determination  of 
Gaussian  wavelet  parameters,  and 

3.1b  (cont’d).  (b)  Misalignment  of  signals  which  produces  sig-  43 
nificant  bias  phase  and  group  delay.  This  arises  because  a 
period  of  25  seconds  is  considered  compared  to  20  seconds 
in  Fig.  3.1a. 

3.2  The  starting  model  (dashed  fine)  used  in  synthetic  test  of  45 
GSDF  inversion  algorithm  and  the 

3.3  The  velocity  time  histories  of  ’observed  seismograms’  46 
(solid  line)  and  the  those  from  the  starting  model  (dashed 
line)  are  both  filtered  in  the  frequency  band  0.01-0.05  Hz 

by  using  a  Butterworth  filter  with  four  poles.  The  plotted 
seismograms  are  normalized  according  the  maximum 
amplitude  of  each  component  in  current  frequency  band. 

It  is  clear  to  see  that  the  starting  model  is  not  close  to  the 
’true’  model  and  this  may  demonstrate  the  ability  of  the 
inversion  programs  to  resolve  structure. 

3.4  After  12  iterations,  the  final  inversion  result  show  an  pre-  47 
dieted  waveforms  almost  identical  to  the  ’observed  seis¬ 
mograms’. 

3.5  The  comparison  between  the  final  model  and  the  ’true  48 
model’.  We  can  see  that  the  2-layer  crust  is  very  close  to 

the  ’true’  model,  but  that  in  the  upper  mantle  the  model 
shows  some  ’zig-zag’  pattern. 
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3.6  The  result  of  inversion  using  the  second  procedure.  This  50 
inversion  procedure  inverts  S  wave  velocity  from  Love 
wave  and  after  that  inverts  P  wave  velocity  from  Rayleigh 
wave  by  assuming  no  anisotropy  effect.  This  shows  a 
acceptable  waveform  match  in  phase,  but  not  in  envelope. 

3.7  The  model  inverted  by  the  second  inversion  procedure.  49 

3.8a  Waveform  fit  of  the  final  inverted  model  for  ALE  in  the  52 
frequency  band  of  0.005-0.03  Hz.  The  signal  which  arrives 
at  1020  seconds  and  around  1300  seconds  are  S  and  SS 
phase,  respectively.  From  the  S  phase  waveform,  we  can 
say  that  the  source  time  function  used  in  this  inversion  is 
a  little  short  but  is  close  enough. 

3.8b  (Cont’d).  Waveform  fit  in  the  0.005-0.05  Hz  frequency  53 
band. 

3.9  The  inverted  models  for  ALE.  Two  models  (PSV  and  SH)  54 
were  inverted  for  Love  wave  and  Rayleigh  Wave.  From 
these  two  models,  there  is  one  slight  anisotropic  zone 
above  200  km. 

3.10  The  final  inverted  model  for  FCC.  There  is  4.5%  55 

anisotropy  effect  exists  between  70  and  140  km. 

3.11a  The  waveform  fitting  for  final  model  in  FCC  are  shown  at  56 
three  frequency  bands  :  (a)  0.01-0.03  (b)  0.01-0.05  (c) 
0.01-0.1  Hz.  The  SS  phase  arrives  at  660  seconds.  The 
inverted  model  can  fit  fundamental  mode  Love  wave  and 
Rayleigh  wave  waveforms  as  high  as  0.1  Hz,  but  it  lacks 
the  ability  to  simulate  the  higher  modes. 

3.11b  (Cont’d).  (b)  0.01-0.05  Hz.  57 

3.11c  (Cont’d).  (c)  0.01-0.1  Hz.  The  SS  phase  arrives  at  660  sec-  58 
onds.  The  inverted  model  can  fit  fundamental  mode  Love 
wave  and  Rayleigh  wave  waveforms  as  high  as  0.1  Hz, 
but  it  lacks  the  ability  to  simulate  the  higher  modes. 

3.12  The  inverted  models  for  PAS.  There  is  no  clear  anisotropy  59 
effects.  The  crustal  Q  is  very  low. 

3.13a  Waveform  fitting  are  shown  on  three  frequency  bands  :  (a)  60 

0.01-0.03  (b)  0.01-0.05  and  (c)  0.01-0.1  Hz.  The  S  wave  sig¬ 
nal  arrives  at  350  seconds. 
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3.13b  (Cont’d).  (b)  0.01-0.05  Hz. 
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3.13c  (Cont’d).  (c)  0.01-0.1  Hz.  62 

4.1  In  GA,  models  are  randomly  generated,  so  there  always  68 
have  some  ’ZIG-ZAG’  patterns.  For  obtaining  a  smooth 
background  velocity  model,  a  smoothing  mechanism  is 
introduced.  Without  changing  the  vertical  shear  wave 
travel  time,  the  smoothed  model  (dashed  line)  has  less 
’ZIG-ZAG’  pattern  than  the  original  model  (solid  line). 

4.2  To  know  how  many  generations  is  necessary  for  surface-  70 
waveform  modeling,  an  experiment  with  large  generation 
nxmiber  (500)  is  tested.  The  result  shows  that  GA  can  find 

a  fairly  good  model  within  50  generations.  After  that,  the 
model  improvement  proceeds  less  rapidly. 

4.3a  Using  GA  to  model  the  surface- waveform  of  station  73 
ANMO.  In  this  test,  two  GA  runs  were  conducted  for  Love 
and  Rayleigh  wave  respectively.  The  searched  models  are 
shown  in  these  figures,  the  heavy  hues  are  the  search 
bounds,  the  thin  lines  eire  searched  model  which  have 
their  goodness-of-fit  greater  than  a  certain  value  shown 
at  the  left  bottom,  and  the  best  model  is  plotted  as  thick 
black  line,  (a)  Searched  models  for  Love  wave. 

4.3b  (cont’d).  (b)  Searched  models  for  Rayleigh  wave.  72 

4.4  For  ANMO,  the  observed  waveforms  (solid  line)  and  pre-  74 
dieted  waveforms  (dashed  line)  which  generated  for  the 
best  searched  models  for  Love  and  Rayleigh  wave  at  fi*e- 
quency  range  0.01  to  0.1  Hz. 

4.5  For  ANMO,  the  filtered  cross-correlations  of  observed  and  75 
synthetic  waveforms  at  different  frequency  bands.  The 
number  at  the  right  of  cross-correlation  traces  is  the 
cross-correlation  value  at  zero-lag  which  is  used  to  con¬ 
struct  the  value  of  goodness-of-fit. 

4.6a  Separated  GA  searchs  for  Love  and  Rayleigh  wave  76 
recorded  at  TUG.  (a)  Searched  models  for  Love  wave. 

4.6b  (cont’d).  (b)  Searched  models  for  Rayleigh  wave.  77 

4.7  A  joint  GA  search  using  both  Love  and  Rayleigh  wave-  78 
forms  at  TUG. 
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4.8a  For  TUC,  filtered  waveforms  at  frequency  range  0.01  to  79 
0.1  Hz  for  observed  and  predicted  seismograms,  (a)  The 
predicted  seismograms  were  generated  using  the  best 
models  from  separated  searched  for  Love  and  Rayleigh 
wave  (Figures  4.6ab). 

4.8b  (cont’d).  (b)  The  predicted  seismograms  were  generated  80 
using  the  best  model  from  a  joint  searched  for  both  Love 
and  Rayleigh  wave  (Figure  4.7). 

4.9a  For  TUC,  the  cross-correlations  for  the  observed  and  S3m-  81 
thetic  seismograms  which  were  generated  using  GA 
searched  models,  (a)  The  synthetics  were  computed  using 
the  best  models  from  separated  searched  models  for  Love 
and  Rayleigh  wave  (Figures  4.6ab). 

4.9b  (cont’d).  (b)  The  synthetics  were  computed  using  the  best  82 
model  from  a  joint  search  for  both  Love  and  Rayleigh 
wave  (Figure  4.7). 

4.10  The  models  fi*om  GA  search  for  Rayleigh  wave  recorded  at  83 
WMOK. 

4.11  The  waveform  comparison  of  observed  and  synthetic  seis-  84 
mograms  for  WMOK  (Figure  4.10). 

4.12  The  cross-correlation  of  the  best  GA  searched  model  for  85 
WMOK  (Figure  4.10). 
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CHAPTER  1 
INTRODUCTION 


1.1  The  Problem 

To  verify  seismological  inferences  about  the  Earth,  matching 
observed  seismograms  an  ultimate  task  for  modern  seismology.  After 
the  success  in  modeling  long-period  seismograms  and  routinely  esti¬ 
mating  earthquake  mechanisms  in  the  early  80s,  seismologists  are 
now  challenged  to  model  entire  regional  seismograms  at  the  high  fre¬ 
quencies.  In  this  study,  two  different  methods  will  be  used  to  model 
the  surface-wave  waveform.  One  of  the  proposed  algorithms  uses  the 
gradient  information  from  the  hypersurface  of  the  misfit  function  to 
formulate  the  inversion  algorithm,  the  other  is  based  on  a  global 
search  technique.  The  advantages  and  weaknesses  of  each  proposed 
algorithm  will  be  examined  in  terms  of  robustness  and  uniqueness  of 
the  resultant  model.  The  ultimate  purpose  of  this  t3q)e  of  research  is 
to  quantify  the  limitations  of  one-dimensional  earth  models  as  an  ini¬ 
tial  step  in  approximating  3-D  structures. 

1.2  Review  of  Related  Literature 

Since  the  birth  of  quantitative  seismology,  seismologists  have 
shared  the  same  goal:  that  one  day  we  can  use  all  geophysical  knowl¬ 
edge  to  design  an  intelligent  system  which  can  routinely  analyze  £dl 
near  real-time  earthquake  data,  predict  future  earthquakes,  and  sys¬ 
tematically  update  our  knowledge  about  the  Earth.  To  reach  this  goal, 
there  are  several  important  tasks  to  accomplish:  collecting  earthquake 


1 
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data,  developing  more  sophisticated  computational  techniques,  and 
finding  a  better  way  to  interpret  data  systematically. 

For  the  past  half  century,  seismological  research  into  Earth  struc¬ 
ture  can  be  roughly  divided  into  three  scales:  planetary  scale  (global 
seismology),  regional  scale  (regional  seismology;  a  few  hundred  kilome¬ 
ters),  and  small  scale  (exploration  seismology;  a  few  tens  of  kilome¬ 
ters).  The  development  for  these  three  scales  has  not  been  uniform. 
Global  and  exploration  seismology  were  developed  more  systematically 
than  regional  seismology.  The  reasons  for  this  phenomenon  are  differ¬ 
ences  in  interests  of  science  and  industry,  in  data  coverages,  and  in  the 
theoretical  computational  methods  in  seismology.  In  this  review,  his¬ 
torical  developments  will  be  briefly  described,  possible  future  develop¬ 
ing  directions  will  be  outlined,  and  the  emphasis  of  this  study  will  be 
addressed. 

1.2.1  Development  Of  Global  Seismology 

Although  earthquakes  happen  in  the  world  everyday,  the  cur¬ 
rently  available  data  are  very  limited  in  both  spatial  and  temporal  dis¬ 
tribution.  Confined  by  the  dataset,  only  the  low-wavenumber  (har¬ 
monic)  global  structimal  features  could  be  estimated.  In  fact,  even 
though  research  in  global  seismology  has  been  boosted  by  specific  sci¬ 
entific  interests,  the  global  earth  structures  determined  is  not  unique. 
However,  global  seismology  developed  very  systematically.  For  exam¬ 
ple,  scientists  inferred  earth  structure  using  different  types  of  data 
(free  oscillations,  surface  waves,  and  body  wave  travel  times)  and  ref¬ 
erence  models  have  been  established  (e.g.  Dziewonski  and  Anderson 
1981).  Based  on  the  reference  models,  routine  analysis  of  seismic 
sources  for  large  earthquakes  became  possible  (Dziewonski,  Chou,  and 
Woodhouse,  1981;  Dziewonski  and  Woodhouse,  1983;  Sipkin,  1986). 
On  the  other  hand,  global  heterogeneity  (obtained  using  tomographic 
techniques)  have  been  well  studied  using  different  data  such  as  phase 
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and  group  velocities  of  long-period  siuface  waves  (Nakanishi  and 
Anderson,  1983,  1984;  Nataf,  Nakanishi,  and  Anderson,  1986),  long- 
period  waveforms  (Woodhouse  and  Dziewonski,  1984),  P-velocities 
from  teleseismic  events  (Dziewonski,  1984),  and  normal  modes  (Mas¬ 
ters  et  al.,  1982).  In  these  studies,  the  horizontal  resolution  is  greater 
than  1000  km  (Zhang  and  Tanimoto,  1992),  and  many  significant  fea- 
tinres,  such  as  mid-ocean  ridges,  require  better  resolution. 

There  are  several  ways  to  increase  the  resolving  power  of  seismo- 
logical  techniques.  The  first  way  is  to  use  more  data  in  inversion. 
Although  the  number  of  global  seismic  stations  are  increasing,  the 
earthquake  dataset  is  increasing  steadily  but  slowly.  A  few  decades 
may  be  needed  to  meet  the  required  size  for  stud3dng  detailed  features. 
Under  such  a  limitation,  seismologists  have  to  find  better  ways  to 
increase  the  usable  data  which  means  to  efficiently  use  the  current 
dataset.  One  possibility  is  to  use  waveforms  data  directly  and  another 
possibility  is  to  use  shorter  period  data  from  smaller  earthquakes.  A 
difficulty  associated  with  using  waveforms  is  the  availability  of  accu¬ 
rate  source  mechanism  of  smaller  regional  earthquakes.  The  second 
way  to  increase  resolution  is  reducing  some  uncertainties  caused  by 
the  inacctmate  shallow  structure.  Mooney  et  al.  (1996)  show  improve¬ 
ment  of  the  resolving  power  in  global  structure  when  using  an 
improved  global  crustal  model  for  station  corrections.  The  third  way 
for  increasing  resolution  is  to  find  new  ways  to  interpret  data.  Two 
experiments  have  been  reported.  One  uses  body  wave  travel-time 
residuals  (Su  and  Dziewonski,  1992),  the  other  uses  differential  travel 
times  (Woodward  and  Masters,  1991).  The  typical  horizontal  resolving 
length  for  such  experiments  is  550  km  (Su  et  al.  1992). 

It  is  clear  that  the  current  effort  on  global  seismology  is  directed 
toward  increasing  resolution.  To  achieve  this,  more  knowledge  is 
required  about  the  regional  structure.  More  complicated  techniques  to 
invert  different  tjqjes  of  data  simultaneously  must  also  be  developed 
and  proven. 
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1.2.2  Development  Of  Exploration  Seismology 

The  purpose  of  exploration  seismology  is  search  for  resources  by 
obtaining  a  high-resolution  image  of  underground  structure.  To  do 
this,  artificial  sources  and  dense  receiver  coverage  are  necessary 
instead  of  relying  on  natural  earthquakes  and  sparse  seismic  stations. 
Huge  exploration  seismic  datasets  enable  seismologists  to  investigate 
the  heterogeneity  of  shallow  structure,  even  only  travel  time  data  are 
used.  In  the  early  stages  of  exploration  seismology,  it  was  natural  to 
use  the  first  arrival  data  (refraction  signal)  to  infer  layer  structures. 
But  the  problem  of  non-uniqueness  of  solutions  has  pushed  the  joint 
use  both  refraction  and  wide-angle  reflection  data.  The  use  of  both 
refraction  and  reflection  data  to  invert  1-D  (e.g.  Braile,  1973),  and  2-D 
structure  (e.g.  Zelt  and  Smith,  1992)  is  well  reported. 

However,  as  pointed  out  by  many  researchers  (e.g.  Huang  et  al. 
1986),  the  precise  identification  of  refraction  and  reflection  signals 
(triplication)  is  crucial  for  the  success  of  such  inversion  schemes.  Dur¬ 
ing  data  processing  steps  like  picking  phases,  artifacts  can  be  intro¬ 
duced.  Sometimes  it  is  worse  when  it  is  necessary  to  interpolate  travel 
time  data  in  order  to  avoid  several  problems  associated  with  data 
sparseness.  All  these  procedures  will  introduce  uncertainties  and 
errors  in  inversion  results.  Another  imperfection  about  this  inversion 
algorithm  is  that  it  can  not  retrieve  the  shear  velocity  information. 
Even  though  a  system  of  signal  enhancement  techniques  such  as  CDP 
stacking  and  migration  were  developed  in  oil  industry,  oil  companies 
still  need  new  techniques  to  improve  the  image  of  the  subsurface  struc¬ 
ture.  All  these  factors  challenge  seismologists  to  find  a  better  way 
solving  problems. 

Waveform  inversion  seems  to  be  a  way  to  avoid  introducing  artifi¬ 
cial  effects  and  to  use  more  information  than  travel  times.  Unfortu¬ 
nately  seismologists  find  it  very  difficult  to  directly  invert  the  wide¬ 
band  seismic  waveform,  and  a  two-step  approach  is  often  used 
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(Nowack  and  Braile,  1993).  Therefore  in  the  exploration  seismology, 
there  are  two  different  problems  regarding  waveform  inversion:  the 
first  is  how  to  retrieve  the  background  velocity  (to  match  the  low- 
frequency  part  waveform),  and  the  other  is  how  to  match  the  entire 
waveform  once  the  background  velocity  is  available.  Much  research  is 
reported  on  fitting  the  entire  waveform  when  a  good  starting  model  is 
available  (e.g.  Mora,  1987;  Pan  and  Phinney,  1989;  and  Crase  et  al. 
1990).  Although  the  assumption  of  availability  of  a  good  starting 
model  seems  suitable  because  of  standard  velocity  analysis  procedures 
widely-used  in  oil  industry,  there  has  been  no  progress  in  retrieving 
background  velocities  using  waveform  inversion  technique  because  of 
its  fully  nonlinear  character. 

Recently,  there  are  two  good  experiments  in  retrieving  back¬ 
ground  velocity  were  performed  by  Keren  et  al.  (1991)  and  Bunks  et  al. 
(1995).  Keren  et  al.  (1991)  try  to  find  the  possible  velocity  model  using 
a  Simulated  Annealing  algorithm.  In  their  experiment,  they  found  the 
existence  of  local  minima  on  the  hypersurface  of  the  misfit  function  is 
caused  by  a  cycle-skipping  effect,  and  a  posteriori  probability  density 
can  be  computed  by  using  simulated  annealing  algorithm.  As  outlined 
previously,  seismologists  know,  from  experience,  it  is  not  a  simple  task 
to  match  the  entire  waveform.  It  seems  that  at  least  two-step 
approach  is  needed.  Based  on  this  idea.  Bunks  et  al.  (1995)  try  to  use 
different  low-pass  filtered  waveform  in  inversion.  They  first  inverted 
for  a  simple  1-D  model  from  the  low-frequency  signal  (0-7  Hz),  then 
inverted  for  a  gross  2-D  model  using  middle-frequency  waveform  (0-10 
Hz),  and  finally  obtained  a  detail  structure  from  inverting  high- 
frequency  waveform  (0-15  Hz).  The  results  of  the  Monte  Carlo  search 
and  multiscale  inversion  are  encouraging.  These  experiments  show 
the  good  resolution  can  we  expect  when  we  have  enough  data  coverage. 
As  commented  by  Nowack  and  Braile  (1993),  the  current  state  of  wave- 
field  inversion  in  exploration  seismology  is  to  find  a  best  inversion  pro¬ 
cedure. 
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1.2.3  The  Possible  Future  Direction  Of  Development 

We  can  see  that  global  seismology  is  currently  trjdng  to  improve 
the  resolution  of  the  tomographic  image  and  determination  of  focal 
mechanism  of  smaller  events  by  extending  Harvard  CMT  method 
(Ekstrom,  1996).  In  exploration  seismology,  high  resolution  images  of 
subsurface  structure  can  be  obtained,  but  only  for  shallow  depths.  So 
there  is  an  existing  gap  between  this  two  scales.  From  the  develop¬ 
ment  of  both  global  seismology  and  exploration  seismology,  we  believe 
the  future  attention  will  focus  on  systematically  stud3dng  the  regional 
structures. 

Theoretically,  we  should  directly  model  the  waveform  for  3-D 
structures,  but  practically,  it  is  impossible  because  of  our  current  lim¬ 
ited  knowledge  about  the  complicated  structures  and  our  sparse  data 
coverage.  So  far,  there  is  no  report  on  successful  modeling  waveform 
for  3-D,  even  2-D  structure.  The  attempt,  such  as  Wdale  and  Helm- 
berger  (1992),  for  matching  the  complicated  waveform  using  2-D  finite- 
difference  method  is  not  successful.  Due  to  the  limited  computing  abil¬ 
ity  of  S3mthetic  seismograms  in  a  realistic  earth,  two  directions  should 
be  considered  for  separated  source  inversion  and  structure  inversion 
regarding  extracting  more  information  in  a  simplified  1-D  model.  The 
first  is  inversion  using  broadband  waveform  data.  The  second  is  inver¬ 
sion  using  as  many  different  types  data  as  possible.  The  reason  for 
this  is  that  while  broadband  waveform  data  place  better  constraints 
than  any  single  data  type,  additional  constraints  can  be  obtained  from 
other  data  such  as  pP  -  P  providing  constraints  on  source  depth.  Once 
seismologists  develop  all  these  tools,  the  final  stage  is  to  build  a  sys¬ 
tem  which  will  automatically  locate  earthquakes,  analyze  source 
mechanisms,  and  indicate  the  "anomalous"  earthquakes  for  further 
study.  All  of  these  will  greatly  contribute  to  other  branch  of  geoscience 
researches  such  as  delineating  the  tectonic  processes,  and  even  moni¬ 
toring  the  failxire  stress  change  in  a  fault  system  (Stein  et  al.  1992).  In 
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the  meantime,  seismologists  should  continuously  develop  forward  and 
inverse  algorithms  for  2-D  and  3-D  structures,  and  a  further  experi¬ 
ment  regarding  the  iterative  source  and  structure  inversion  could  be 
tested. 

Recently,  there  has  been  much  research  in  focal  mechanism  analy¬ 
sis  using  regional  broadband  waveform.  This  work  using  different 
subsets  data  such  as  surface-wave  (Thio  and  Kanamori,  1995),  body- 
wave  waveform  (Dreger  and  Helmberger,  1993),  full  waveform  (Mao, 
Panza,  and  Suhadolc,  1994),  and  even  one  single  station  (Fan  and  Wal¬ 
lace,  1991).  Zhou  et  al.  (1994)  and  Zhu  et  al.  (1996)  showed  one  inter¬ 
esting  method  which  uses  both  high-frequency  body-wave  waveform 
and  low-ffequency  full  waveform  in  searching  for  the  best  focal  mecha¬ 
nism.  In  this  study,  we  will  concentrate  on  the  regional  structure 
inversion  problem  since  earth  structure  is  fundamental  to  source 
parameter  estimation. 

1.2.4  Emphasis  Of  This  Study 

To  model  regional  broadband  waveforms,  there  are  two  different 
approaches.  The  first  approach  is  to  match  different  phases  of  body- 
wave  such  as  P„;  and  S^i-  Zhao  and  Helmberger  (1991)  proposed  one 
strategy  to  model  the  waveform.  First,  they  model  long-period  data 
(P„  and  Sn)  assuming  a  single  crustal  layer  model  and  get  the  average 
crustal  and  upper  mantle  P,  S  velocities,  and  crustal  thickness.  Next, 
they  model  the  fundamental  mode  Rayleigh  wave  by  adding  some  lay¬ 
ers  in  the  crust  to  the  previous  simple  model.  Third,  they  model  those 
short-period  waveforms  riding  on  P„  and  by  adding  some  layers  in 
upper  mantle  structure.  Finally,  they  do  the  final  adjustments  by  trial 
and  error.  This  is  a  great  strategy  which  utilizes  many  body-wave 
phases  with  clear  physical  meanings,  but  it  requires  an  experienced 
seismologist.  The  second  approach  requires  matching  the  surface- 
wave  waveform.  One  example  for  this  approach  is  Hnear  inversion 
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method  shown  by  Gomberg  and  Masters  (1988).  They  use  modal 
superposition  to  create  synthetic  seismograms  and  perturb  them  with 
respect  to  layer  parameters  (such  as  S-velocity).  Using  these  partial 
derivatives,  they  formulate  a  linear  inverse  problem  for  the  differen¬ 
tial  seismogram  between  observed  and  synthetic.  In  this  study,  we 
will  develop  other  inversion  methods  for  sinface-waveform  because 
this  approach  has  the  potential  to  be  a  fully  automatic  process.  In 
most  regions,  the  sparse  data  coverage  does  not  support  directly  mod¬ 
eling  waveform  in  2-D,  or  3-D  structures.  Therefore,  the  S3nithetic 
seismograms  used  in  this  study  will  be  generated  for  a  1-D  model 
using  mode  summation  method  (Wang,  1981). 

Using  a  cross-correlation  technique,  Lerner-Lam  and  Jordan 
(1983)  formulated  an  inversion  algorithm  to  obtain  earth  structure 
from  long-period  teleseismic  surface- wave  waveforms.  A  successful 
application  was  presented  by  the  same  group  (Lerner-Lam  and  Jordan 
1987)  on  investigating  the  continental  thickness  of  Eurasia.  Similar  to 
this  idea,  Cara  and  Leveque  (1987)  showed  a  slightly  different  way  to 
extract  the  "secondary  observations"  and  formulated  an  algorithm  to 
invert  these  "secondary  observations".  They  used  teleseismic  seismo¬ 
grams  which  have  well  separated  fundamental  mode  and  higher 
modes  to  ensure  their  success.  However,  seismograms  show  more  com¬ 
plicated  mode  interference  effects  on  the  regional  scale. 

Gee  and  Jordan  (1992)  introduced  Generalized  Seismological  Data 
Functionals  theory  (GSDF)  to  deal  with  surface-wave  modal  interfer¬ 
ence  effects.  In  this  study,  an  algorithm  based  on  the  GSDF  theory  is 
proposed  and  will  be  tested  using  an  earthquake  that  occurred  in  the 
west  Texas  and  was  well  recorded  throughout  North  America.  The 
comparison  between  GSDF  and  linear  inversion  (Gomberg  and  Mas¬ 
ters  1988)  will  also  be  addressed. 

Using  linear  inversion  theory,  a  good  starting  model  is  necessary 
to  guarantee  the  success  of  final  convergence.  It  is  necessary  to  use  a 
global  search  method  to  ease  the  dependence  on  apriori  information. 
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A  second  algorithm  is  considered.  To  use  Genetic  Algorithms  (GA)  in 
searching  models  which  produce  good  surface-wave  waveform  fit. 
Hopefully,  these  two  algorithms  will  provide  more  insight  in  how  to 
build  a  more  intelligent  and  robust  inversion  algorithm. 


CHAPTER  2 

DATA  SET  USED  FOR  COMPARISON 


In  this  study,  the  different  algorithms  will  be  tested  using  the 

April  14,  1995  western  Texas  earthquake  data.  This  earthquake  is 

well  recorded  by  broadband  seismometers  in  the  northern  America. 

The  currently  data  set  is  collected  from  IRIS  (Incorporated  Research 

•» 

Institutions  for  Seismology),  CNSN  (Canadian  National  Seismological 
Network),  USGS  (United  States  Geological  Survey),  UNAM,  and 
PASSCAL  experiments.  The  station  distribution  map  is  shown  as  Fig- 
ure2.1,  and  the  station  information  is  listed  on  Table  2.1.  In  Table  2.1, 
the  information  of  each  station’s  location,  azimuth,  hypocenter  dis¬ 
tance,  P  wave  first  arrival  polarity,  takeoff  angle,  and  whether  the  dis¬ 
persion  data  used  in  searching  Rayleigh  or  Love  wave  radiation  pat¬ 
tern  are  shown. 

Using  these  data,  the  different  algorithms  for  modeling  surface- 
waveform  will  be  tested.  Before  this  stage,  the  source  parameters  will 
be  reestimated,  i.e.  source  depth  and  focal  mechanism.  Table  2.2  lists 
three  different  source  parameters  which  include  Harvard  CMT  solu¬ 
tion,  USGS  Sipkin’s  solution,  and  the  reestimated  source  parameters 
used  in  this  study. 

Comparing  USGS  and  Heirvard  CMT  solution,  they  both  have 
good  agreement  in  h5q)ocenter  location  and  origin  time.  However,  the 
focal  mechanism  is  very  different  from  these  two  solutions.  From 
experiences,  the  published  source  depth  for  shallow  events  usually 
have  huge  uncertainty.  Therefore,  we  simply  use  USGS’s  h3q)ocenter 
location  and  origin  time  and  will  reestimate  the  focal  mechanism  and 
depth. 
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Table  2.1  The  station  informations 


STA 

LON 

LAT 

AZ 

DIST 

P 

R 

L 

IH-: 

ALE 

-62.3500 

82.5033 

6.0800 

6036.3398 

-1 

0 

0 

29.1 

ANMO 

-106.4567 

34.9462 

331.3424 

596.7890 

-3 

1 

1 

71.1 

ARC 

-124.0750 

40.8770 

307.8156 

2210.4399 

0 

1 

0 

47.1 

BAR 

-116.6720 

32.6800 

285.3669 

1295.2960 

0 

0 

1 

65.1 

BBB 

-128.1133 

52.1847 

327.2636 

3169.0320 

0 

1 

1 

37.1 

BINY 

-75.9861 

42.1993 

53.9871 

2774.6509 

0 

1 

1 

39.1 

BKS 

-122.2350 

37.8770 

300.9237 

1934.0950 

-3 

1 

0 

57.. 

BMN 

-117.2218 

40.4315 

315.6311 

1689.8090 

-1 

1 

1 

60.1 

CAIG 

-100.2680 

17.0480 

167.3191 

1496.1100 

0 

1 

1 

62.: 

CALB 

-118.6270 

34.1430 

290.0596 

1503.6780 

-5 

1 

0 

62. 

CBKS 

-99.7374 

38.8140 

18.1580 

1004.1560 

0 

1 

0 

mSSJl 

CCM 

-91.2446 

38.0557 

48.9178 

1408.1630 

0 

1 

1 

63.: 

CEH 

-79.0928 

35.8908 

68.1402 

2340.4871 

-3 

0 

0 

43. 

CMB 

-120.3850 

38.0350 

303.3337 

1789.8700 

-1 

1 

1 

60.i 

COL 

-147.7933 

64.9000 

334.6238 

4911.2661 

-1 

0 

0 

32. 

COR 

-123.3032 

45.5857 

317.7288 

2365.8679 

0 

1 

1 

42. 

CUIG 

-99.1780 

19.3290 

159.9880 

1281.2321 

0 

1 

1 

65. 

CYF 

-109.8660 

37.5542 

324.9922 

1009.3980 

-1 

1 

1 

69. 

DAWY 

-139.4320 

64.0530 

336.5904 

4525.5771 

0 

0 

0 

33. 

DGR 

-117.0090 

33.6500 

289.7421 

1345.7460 

3 

1 

1 

64. 

DLBC 

-130.0272 

58.4372 

334.7115 

3733.0969 

0 

1 

1 

■iftl 

DRLN 

-57.5042 

49.2560 

47.8932 

4368.4370 

0 

0 

0 

34. 

DUG 

-112.8133 

40.1950 

324.5603 

1398.2980 

0 

1 

1 

EDM 

-113.3500 

53.2220 

345.1403 

2677.0239 

-1 

1 

1 

EYMN 

-91.4950 

47.9460 

23.9364 

2208.1650 

-1 

1 

1 

47. 

FCC 

-94.0870 

58.7620 

9.8789 

3244.3350 

-1 

1 

1 

36. 

FFC 

-101.9783 

54.7350 

1.8890 

2719.7180 

-1 

0 

0 

40. 

FRB 

-68.5470 

63.7470 

23.2040 

4453.6089 

-1 

0 

0 

33. 

GAC 

-75.4783 

45.7033 

46.8810 

2962.3689 

0 

1 

0 

38. 

GAR 

-114.1020 

38.8808 

316.9895 

1373.5229 

-1 

1 

1 

63. 

GOL 

-105.3711 

39.7003 

350.4675 

1063.6190 

-5 

0 

0 

mm 

GSC 

-116.8083 

35.3028 

297.4232 

1397.5800 

0 

1 

1 

63. 

HOPS 

-123.0723 

38.9935 

303.5273 

2047.1680 

0 

0 

0 

53. 

HRV 

-71.5580 

42.5060 

55.5348 

3132.9500 

0 

1 

1 

37. 

INK 

-133.5200 

68.3070 

343.7640 

4662.9888 

-1 

0 

0 

33. 

ISA 

-118.4733 

35.6633 

296.9282 

1535.3850 

0 

0 

0 

61. 

JFWS 

-90.2488 

42.9149 

36.0724 

1823.3640 

0 

1 

1 

59. 

JRSC 

-122.2376 

37.4038 

299.4258 

1916.7350 

0 

0 

0 

57. 

LBNH 

-71.9259 

44.2401 

51.8728 

3160.6780 

0 

1 

1 

37. 

LMN 

-64.8060 

45.8520 

51.6319 

3749.1060 

0 

1 

1 

■na 

LMP 

-111.1960 

37.9783 

321.8398 

1121.4580 

-1 

1 

1 

67. 

LMQ 

-70.3267 

47.5483 

46.3837 

3406.0991  . 

0 

1 

1 

36. 

LSCT 

-73.2244 

41.6784 

56.5500 

2975.6770 

0 

1 

1 

38. 

MBC 

-119.3600 

76.2420 

354.7885 

5187.3291 

0 

0 

0 

32. 

MDW 

-112.3570 

38.8927 

321.6075 

1264.8280 

-3 

1 

1 

65. 

MHC 

-121.6420 

37.3420 

299.7602 

1864.6530 

-1 

0 

0 

58. 

MIAR 

-93.5730 

34.5457 

60.0795 

1032.6490 

5 

1 

0 

68. 

MIN 

-121.6050 

40.3450 

308.9899 

1997.3640 

-1 

1 

0 

55. 

MLAC 

-118.8340 

37.6310 

303.8644 

1647.0140 

-3 

1 

1 

mm 

MLC 

-116.7680 

38.9883 

311.8189 

1564.3490 

-1 

1 

1 

61. 

MNV 

-118.1530 

38.4328 

307.6257 

1634.5140 

-1 

0 

0 

61. 

NEE 

-114.5960 

34.8230 

298.4929 

1172.2111 

0 

1 

1 

67. 

NEW 

-117.1200 

48.2630 

333.4136 

2318.1299 

-1 

1 

1 

44. 

NWC 

-113.5566 

37.6332 

313.6112 

1248.7170 

0 

0 

0 

66. 

ORV 

-121.5000 

39.5560 

306.8155 

1949.1410 

-1 

1 

1 

57. 

OXIG 

-96.7230 

17.0820 

153.8263 

1606.6180 

0 

0 

0 

61. 

PAS 

-118.1717 

34.1483 

190.9800 

1462.8500 

0 

1 

1 

62. 

b5i>OOObOWOCnCOC;iOO;OOOOOlO>COOlh-‘CX)0500rf^COCO<ia>OOOiOOM<IOiCOt^005CO<IOaiOOOO<ICOOH*tOO>»I^^Ot0050Cn 
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PFO 

-116.4553 

33.6092 

290.0471 

1294.6429 

0 

0 

0 

65.2 

PGC 

-123.4508 

48.6500 

325.7721 

2660.2791 

0 

1 

1 

40.4 

PH2 

-114.9839 

37.6645 

310.4570 

1352.8110 

-3 

1 

1 

64.2 

PKDl 

-120.4249 

35.8892 

295.8760 

1711.1760 

-1 

1 

0 

60.5 

PLIG 

-99.5010 

18.3890 

162.7272 

1370.8290 

1 

1 

1 

63.9 

PMB 

-123.0765 

50.5188 

329.3981 

2786.8779 

-1 

1 

1 

39.8 

PNIG 

-98.1270 

16.3950 

159.8122 

1624.5060 

0 

1 

1 

61.0 

PNT 

-119.6200 

49.3200 

331.5872 

2522.8210 

0 

1 

0 

40.9 

RCC 

-110.5850 

40.5163 

331.9426 

1313.9041 

-1 

1 

1 

64.8 

RES 

-94.9000 

74.6870 

317.4068 

4965.4668 

-1 

0 

0 

32.6 

RPV 

-118.4035 

33.7438 

289.0763 

1474.4340 

-3 

1 

1 

62.4 

RTS 

-115.7900 

39.6681 

315.9504 

1541.3740 

-3 

1 

0 

61.7 

SADO 

-79.1417 

44.7694 

46.0372 

2659.1970 

-1 

1 

1 

40.4 

SAG 

-121.4450 

36.7650 

298.0289 

1827.4780 

-3 

0 

0 

59.5 

SCHQ 

-66.8336 

54.8319 

36.0646 

3975.2800 

-3 

1 

1 

34.7 

SMTC 

-115.7200 

32.9490 

287.4082 

1212.3929 

0 

1 

0 

66.7 

SRS 

-110.6010 

38.9126 

327.1899 

1168.0090 

-1 

1 

1 

67.3 

SSPA 

-77.8914 

40.6401 

56.4481 

2567.4299 

-1 

0 

0 

40.6 

TUG 

-110.7846 

32.3096 

289.6425 

745.3420 

-2 

1 

1 

69.6 

ULM 

-95.8750 

50.2498 

13.5989 

2305.5010 

-2 

1 

1 

44.6 

use 

-118.2870 

34.0210 

290.3445 

1470.1250 

0 

1 

1 

62.5 

VTV 

-117.3330 

34.5670 

293.5786 

1399.8250 

0 

1 

0 

63.4 

WALA 

-113.9115 

49.0586 

339.7463 

2270.9600 

-1 

1 

1 

45.5 

WCP 

-114.1670 

40.5242 

322.1671 

1502.9880 

-1 

1 

1 

62.1 

WDC 

-122.5400 

40.5800 

308.5890 

2079.6699 

-3 

1 

1 

52.7 

WHY 

-134.9906 

60.6597 

334.5919 

4102.5229 

-1 

1 

1 

34.4 

WMOK 

-98.7810 

34.7379 

39.5295 

655.0750 

3 

1 

0 

70.3 

WMT 

-111.8380 

40.1111 

327.0212 

1338.3270 

0 

1 

1 

64.4 

WVOR 

-118.6367 

42.4339 

318.8431 

1921.1140 

-1 

0 

0 

57.7 

YBH 

-122.7195 

41.7318 

311.4840 

2154.2839 

0 

1 

1 

49.8 

YKW3 

-114.6050 

62.5620 

350.4589 

3681.6750 

-1 

1 

0 

35.9 

ZIIG 

-101.4650 

17.6070 

171.9038 

1414.2450 

0 

0 

0 

63.2 

From  stations  which  have  source-receiver  distance  greater  than 
3600  km,  it  is  very  clear  to  see  a  general  7-9  seconds  time  difference 
between  P  and  pP  phases  (Figure2.2).  At  such  distance,  the  observed 
first  arrival  P  is  beyond  the  upper  mantle  P  triplication  zone  which 
caused  by  upper  mantle  discontinuities.  This  means  that  source  depth 
is  not  as  shallow  as  suggested  by  Harvard  and  USGS  solutions.  For  a 
simple  crust  model,  such  a  pP-P  difference  will  correspond  to  an 
earthquake  deeper  than  20  kilometers. 

The  source  mechanism  used  in  this  study  is  based  on  the  grid 
search  for  Love  and  Rayleigh  wave  radiation  pattern  (Herrmann, 
1974)  using  observed  spectrum  amplitude  of  fundamental  mode  group 
velocities  at  different  frequencies.  Based  on  the  grid  search  for  focal 
mechanism  at  different  source  depths,  the  optimum  result  is  a  point 
source  at  23  km  in  depth,  with  strike,  dip,  rake  as  (114°,  64°,  -101°), 
and  with  seismic  moment  3.0E24  dyne-cm. 
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Fig.  2.1.  Station  distribution  map.  The  stations  are  shown  as  triangles. 
The  focal  mechanism  used  in  this  study  is  shown  as  the  beach  ball. 
This  map  is  generated  using  GMT  (Wessel  and  Smith,  1991) 

From  the  station  distribution  map  (figure  2.1),  we  can  see  that 
most  of  broadband  stations  are  in  the  western  United  States.  The 
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Table  2.2  Source  parameters  for  different  solutions 


Lon 

Lat 

Depth 

Origin  Time 

Strike 

Dip 

Rake 

moment 

HARVAKD 

-103.32 

30.24 

15.0 

00:32:54,2 

117 

53 

-87 

3.82E24 

USGS 

-103.327 

30.261 

13.0 

00:32:55.04 

136 

60 

-86 

4.00E24 

this  study 

-103.327 

30.261 

23.0 

00:32:55.04 

114 

64 

-101 

3.00E24 

Fig.  2.2.  For  teleseismic  records,  a  general  7-9  seconds  travel  time  difference  for  pP-P 
is  observable. 

second  best  cover  region  is  the  northeast  United  States.  The  funda¬ 
mental  mode  Rayleigh  and  Love  wave  spectrum  amplitudes  used  in 
radiation  pattern  search  also  have  more  dense  data  points  in  these  two 
regions.  Looking  at  the  observed  spectrum  amplitude  data,  it  delin¬ 
eates  a  clear  nodal  plane  for  Love  wave  radiation  pattern.  All  three 
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Fig.  2.3.  The  Love  wave  radiation  patterns  for  the  preferred  focal  mechanism  at  23 
km  depth.  The  bars  indicate  attenuation  corrected  spectral  amplitudes  in  cm-sec  nor¬ 
malized  for  geometrical  spreading  to  1000  km. 


solutions  can  provide  good  Love  wave  radiation  pattern  to  match  the 
observed  data.  But  from  Rayleigh  wave  radiation  pattern  fit,  the  rees¬ 
timated  focal  mechanism  is  better  than  other  solutions.  For  compari¬ 
son,  figure  2.3,  2.4,  and  2.5  will  show  the  Love  wave  radiation  pattern 
at  nine  periods  for  the  mechanism  used  in  this  study.  Harvard  CMT 
solution,  and  USGS  Sipkin’s  solution.  Similar  to  Love  wave,  Figure 
2.6,  2.7,  and  2.8  will  display  the  Rayleigh  wave  radiation  patterns. 
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Fig.  2.4.  The  Love  wave  radiation  patterns  for  Harvard  CMT  solution  at  15  km 
depth. 


P-wave  polarities  provide  another  property  to  check  the  focal 
mechanisms.  The  P-wave  polarity  reading  is  shown  in  table  2.1.  The 
negative  number  indicates  the  downward  movement,  and  the  positive 
number  indicates  the  upward  movement.  The  number  1  indicates  a 
impulsive  P-wave  arrival;  the  larger  numbers  indicate  the  emergent 
arriveds  with  less  confident.  The  P-wave  polarities  for  these  three 
focal  mechanisms  are  shown  as  figure  2.9,  2.10,  and  2.11.  In  these 


Fig.  2.5.  The  Love  wave  radiation  patterns  for  Sipkin’s  USGS  solution  at  13  km 
depth. 

figures,  the  cycles  represent  the  impulsive  upward  P  arrivals  while  the 
plus  signs  represent  the  emergent  upward  P  arrivals;  in  the  similar 
way,  the  triangle  represent  the  impulsive  downward  P  arrivals  and  the 
minus  signs  are  for  the  emergent  downward  P  arrivals.  The  reesti¬ 
mated  focal  mechanism  which  have  only  one  inconsistent  station  is 
better  than  the  other  two  solutions. 
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Fig.  2.6.  The  Rayleigh  wave  radiation  patterns  for  the  preferred  focal  mechanism  at 
23  km  depth. 

From  all  these  comparisons,  the  depth  phases,  Rayleigh  and  Love 
wave  radiation  patterns,  and  the  first  arrival  P  wave  polarities,  we  are 
sure  that  the  reestimated  focal  mechanism  is  better  than  the  others. 


Fig.  2.7.  The  Rayleigh  wave  radiation  patterns  for  Harvard  CMT  solution  at  15  km 


F 
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Fig.  2.9.  This  figure  shows  the  observed  P-wave  polarities  and  the  reestimated  focal 
mechanism  which  will  be  used  for  testing  several  surface-waveform  modeling  algo¬ 
rithms. 
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CHAPTER  3 

REVIEW  OF  GSDF  THEORY  AND  FORMULATE 
STRUCTURE  INVERSION  ALGORITHM  USING  GSDF 


3.1  Introduction 

After  successful  observing  the  apparent  anisotropy  of  Eurasian 
upper  mantle  (Gee  and  Jordan,  1988),  Gee  and  Jordan  (1992)  intro¬ 
duced  the  theory  of  Generalized  Seismological  Data  Functionals 
(GSDF)  which  deals  with  the  interference  effects  of  different  surface 
wave  modes  at  a  certain  frequency.  This  theory  allow  seismologists  to 
precisely  measure  the  group  and  phase  delays  between  observed  and 
synthetic  seismograms  for  different  frequencies  even  when  the  modal 
interference  effect  can  not  be  ignored.  The  applications  have  been 
reported  by  Gaherty  et  al.  (1995,  1996)  who  use  this  method  to  invert 
the  upper  mantle  velocity  structure  and  anisotropy.  In  this  study,  we 
extend  the  Generalized  Seismological  Data  Functionals  (GSDF)  theory 
to  the  inversion  of  broadband  waveforms  modeled  as  a  superposition  of 
surface-wave  modes  (Wang,  1981).  We  demonstrate  the  utility  of  this 
inversion  algorithm  by  conducting  a  simple  synthetic  test,  and  then 
apply  the  algorithm  to  real  data  to  see  what  new  knowledge  about 
earth  structure  can  be  obtained  from  this  new  inversion  algorithm. 

3.2  Theory  of  Generalized  Seismological  Data  Func¬ 
tionals 

To  use  the  GSDF  approach,  we  need  to  construct  a  synthetic  seis¬ 
mogram  (s),  an  isolation  function  (/*)  and  single-mode  seismograms. 
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An  isolation  function  is  a  sum  of  single-mode  seismograms  which  may 
represent  the  dominant  part  of  the  observed  seismogram  (s).  Using  a 
cross-correlation  technique,  we  can  quantify  the  similarity  observed 
and  synthetic  seismograms.  For  a  model  which  does  not  significantly 
deviate  from  actual  earth  structure,  the  peak  of  cross-correlagrams 
will  be  located  near  zero  lag-time  for  all  frequency  bands  and  for  all 
windowed  segments  of  the  seismograms.  The  cross-correlagrams 
reflect  the  degree  that  the  model  used  to  create  synthetic  seismogram, 
isolation  function  and  single-mode  seismograms  is  different  from  real 
earth  structure.  To  utilize  this  information,  we  have  to  extract  infor¬ 
mation  from  different  frequency  ranges  and  interpret  it  in  term  of  dif¬ 
ferences  in  earth  structure. 

Gee  and  Jordan  (1992)  call  the  extracted  information  Generalized 
Seismological  Data  Functionals  (GSDF).  In  the  following,  we  give  a 
brief  review  of  how  the  GSDF  is  defined,  how  the  GSDF  is  related  to 
familiar  physical  quantities,  and  how  we  extend  this  theory  to  invert 
for  Earth  structme.  Equations  directly  adapted  from  Gee  and  Jordan 
(1992)  will  be  cited  as  (GJ.#),  where  #  indicates  the  original  equation 
number  in  their  paper. 

First,  we  review  the  definition  and  computation  of  a  GSDF.  We 
can  use  a  Gaussian  wavelet  to  approximate  the  filtered  and  windowed 
cross-correlagrams  at  a  specified  frequency.  Specifically,  the  win¬ 
dowed,  filtered  cross-correlagram  of  an  isolation  function  with  a  seis¬ 
mogram  (either  observed  or  synthetic)  is  modeled  using  a  five- 
parameter  Gaussian  wavelet: 

FiMVCfsit)  »  git)  =  AGa[cr,(^  -  f^)]  cosK(<  -  tp)]  (GJ.5) 

FiWCfsit) «  git)  =  AGaldsit  -  tg)]  cos[&sit  -  tp)]  (GJ.8) 

F^WC^(^)  =  GalcTft]  cos[d5^^] 


where 
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Fj  is  a  Gaussian  frequency  content  filter  with  center  fi*equency  coi 
and  half-bandwidth  (jj, 

W  =  Ga[(T^(^  -  +  0. 01[<7u,(^  -  - 1  is  a  temporal  window 

with  half-bandwidth  (rjj,and  centered  at  usually  =  t^,  where 
tc  is  the  lag-time  of  the  peak  of  FjWC^, 

Cfs{t)  =  fit)  <a  s(t),  0  represents  cross-correlation,  is  the  cross¬ 
correlation  with  the  observed  time  history, 

Cfs(t)  =  fit)  ^  sit),  is  the  cross-correlation  with  the  synthetic  time 
history, 

Ga[x]  s  exp(-  -— ),  Ga  is  the  Gaussian  function, 

JU 

A  is  the  amplitude  of  Gaussian  envelope, 

(7  is  the  half-bandwidth  of  envelope  spectrum, 

CO  is  the  angular  frequency  of  oscillating  wavelet, 
tg  is  the  envelope  group  delay  from  zero  lag-time, 
tp  is  the  wavelet  phase  delay  from  zero  lag-time, 

the  subscript  s  denotes  observed  seismogram,  subscript  s  com¬ 
bined  with  ~  denotes  synthetic  seismogram,  and  the  subscript  f 
with  ~  denotes  isolation  filter. 

From  (GJ.5)  and  (GJ.8),  Gee  and  Jordan  define  four  data  functionals  to 
characterize  the  agreement  between  observed  and  predicted  time  his¬ 
tories: 


1 

II 

(GJ.9) 

^g 

1 

(GJ.IO) 

5t,= - [InA-lnA] 

(GJ.ll) 

Sta=-^[cOs-&f\ 

(GJ.12) 
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These  four  data  functionals  are  related  to  differential  phase  delay,  dif¬ 
ferential  group  delay,  differences  in  logarithmic  amplitudes  and  the 
differences  in  center  frequencies.  These  four  data  functionals  are 
defined  from  two  filtered,  windowed  cross-correlagrams,  they  can  be 
transformed  into  the  more  familiar  quantities  such  as  phase  velocity, 
group  velocity,  and  attenuation  corrections. 

We  will  briefly  describe  this  transformation.  We  note  that  all  fil¬ 
tered  and  windowed  cross-correlagrams  will  be  normalized  by  scaling 
F^WC^  to  unit  amplitude.  Also,  if  the  observed  and  synthetic  seismo¬ 
grams  are  composed  only  of  a  single  mode,  then  the  GSDF  are  easily 
interpreted  in  terms  of  differences  in  modeled  and  actual  phase  veloc¬ 
ity,  group  velocity  and  Q.  For  multi-mode  time  histories  there  inter¬ 
pretation  is  much  more  difficult. 

Before  relating  GSDF  to  physical  quantities,  we  must  know  the 
roles  these  quantities  play  in  wave  propagation.  For  a  known  instru¬ 
ment  response  and  source,  the  difference  between  the  spectrum  of  iso¬ 
lation  filter  f(o))  -  I(co)P(co)S((o)  and  its  corresponding  component  of 
the  observed  seismograms  fico)  =  I((o)P(ct))S(co)  is  called  "differential 
propagation".  They  are  related  by  the  equation 


fico)  =  fim) 
Pico) 


=  Dico)f{(o) 


(GJ.42) 


_  ^iSk(.m)R =  giik(o))-k(a>)]R 


where  R  is  source-receiver  distance.  To  a  first  order  approximation, 
the  differential  propagation  can  be  approximated  as 

r  \ 


Did))  ~  exp 


^0)jSTgio}j)  -ico-  (Oj)STaicOj)j 


(GJ.44) 


-t-exp 


COjSTpicOj)  +  (ct}-  a)j)SZqicOj) 


■N 

y 


The  propagation  effects  of  the  isolation  filter  fit)  and  the 
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corresponding  waveform  fit)  on  the  observed  seismogram  are  defined 
in  four  equations  which  are  described  in  terms  of  familiar  physical 
quantities  such  as  total  phase  delay,  total  group  delay,  and  attenuation 
factor. 


I 

s 

111 

(GJ.45) 

STgiWj)  =  tgicOj)  -  fgicOj) 

(GJ.46) 

StgicOj)  =  ^  fpicoj)[Q~^  -  Q  "^1 

(GJ.47) 

Sr  aim  j)  =  i  Tgi(Oj)[Q~^  -  Q  "^] 

(GJ.48) 

where 

tpicoj)  is  total  phase  delay  of  fit)  at  an  arbitrary  frequency  coj 
with  respect  to  the  original  time, 

Tpicoj)  is  total  phase  delay  of  fit)  at  coj, 

Tgicoj)  is  total  group  delay  of  fit)  at  coj, 
fgicoj)  is  total  group  delay  of  fit)  at  coj, 

Q  is  attenuation  factor  of  fit), 

Q  is  attenuation  factor  of  fit). 

This  means  that  whenever  we  are  able  to  measure  four  differential 
quantities  from  observed  and  synthetic  seismograms  at  a  certain 
frequency  coj,  we  can  approximate  the  true  waveform  spectrum  behav¬ 
ior  fioj)  at  any  arbitrary  frequency  co  to  a  first  order  precision. 

The  GSDF  St,^  based  on  the  filtered,  windowed  cross-correlation,  and 
the  Stx  is  related  to  the  time  domain  windowing  function.  Their  rela¬ 
tionships  are  given  by 

Stg  ~  STgi&f)  +  (1  -  ^\){tc  -  Sxgi&f)]  (GJ.56) 

Stp  ~  STpi&f)  +  (1  -  ^i)(  — ^  —  )[^e  -  dZgi&f)] 


(GJ.57) 
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6ta~  ^ISTaimf)  (GJ.58) 

at,  =  5r,[S,;)  -  (1  -  (GJ.69) 

where 

Sc  is  the  frequency  from  C^, 
is  a  window  width  factor, 

tc  is  the  lag-time  of  the  peak  of  cross-correlagram  FjWC^. 

There  is  an  additional  assumption  that  must  be  stated.  Due  to  the  dif¬ 
ficulty  in  relating  the  isolation  filter’s  corresponding  feature  in 
observed  seismogram,  it  is  not  easy  to  evaluate  F^WC^.  However,  if 
the  windowing  procedure  is  effective  in  isolating  f(t)  from  other 
phases  on  the  seismograms,  then  FjWC^s  »  FjWC^. 

We  can  estimate  tc  from  an  isolated  waveform.  For  non-isolated 
waveforms,  the  Cff  can  not  be  simply  replaced  by  and  this  is  more 
general  situation.  For  non-isolated  waveforms,  the  Gaussian  wavelet 
can  be  represented  as  a  sum  of  interference  effects  of  different  modes. 

F,WC^,(^)  =  i  X  FiWC,„„(^)  (GJ.63) 

m-1  n=l 


FjWC^(«)  =  I  i  F.WC^tf) 

m=l  n=l 

where 

C^n(f)  =  [mode  m  of  fit)]  [mode  n  of  fit)], 

C,nn(t)  =  [mode  m  of  fit)]  is  [  mode  n  of  fit)]. 

We  now  characterize  these  by  the  relation 

FiWC,„„(0  »  A„,„Ga[<?;.(^ -  f^^)]  cos[Smn(.t - 1^")]  (GJ.68) 

where 

Amn=expi-Sft^’'), 
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^  ^  ^ 

(Omn=^f  ~  a  ’ 

In  this  case,  we  need  to  know  four  time  shifts  that  describe  the  devia¬ 
tions  of  FjWCys  from  FjWC^.  Assume  that  the  windowing  and  filter¬ 
ing  effectively  suppress  the  bandwidth  variations,  so  that  dg-  df.  The 
four  time  shifts  are  the  phase  delay  t  p,  a  group  delay  t  g,  and  two 
amplitude  parameters 

t~=-  —  \nA  (GJ.60) 

Of 

By  defining  the  following  notation, 

- tg)]  (GJ.73) 

z;ri/2 

<l>mn  =  («/  -  c}tT)itT  -  -  ia)itp  -  tg)  (GJ.74) 

the  perturbation  expansions  can  be  simplified  by  using  a  set  NxN 
matrices 


(G)„in  -  -^mra  COS 

(GJ.75) 

(®)/nn  ~  Bffiji  sin  (pffiji 

(GJ.76) 

(GJ.77) 

(Sx)mn  =  -  t^)sin(l)r^ri 

(GJ.78) 

Gee  and  Jordan  (1992)  gave  the  following  linearized  relationships 
between  the  GSDFs  and  computable  quantities  from  individual  mode 

branches. 

(5^,,  =  1  •  C  •  -1- 1  •  S  •  ^tp 

(GJ.84) 

=  -  1  •  S  •  ^tq -1-  1  •  C  •  Jtp 

(GJ.85) 

dta=-l-  (Ca  +  Sg)  •  ^tq  -h  1  •  (Cg  -  Sa)  •  Jtp 

(GJ.86) 
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+  1  •  C  •  +  1  •  S  •  ^tg 


Stg=-1-  (Cg  -  Sa)  (Ca  +  Sg)  •  ^tp  (GJ.87) 


-  1  •  S  •  (Jta  +  1  •  C  ■  <Jtg 

The  1  is  a  A^-dimension  vector,  each  element  of  which  equals  one.  The 
St^  is  a  AT-vector  whose  n’th  component  is  St^,  the  perturbation  corre¬ 
sponding  to  the  quantity  x  of  n’th  mode.  In  the  following  section,  we 
will  translate  this  ^t^  into  which  directly  relate  to  seismogram  and 
has  clear  physical  meaning. 


3.3  Structure  Inversion 

To  apply  GSDF  theory  in  structure  inversion,  we  must  construct 
the  inversion  kernel,  G  for  our  inversion.  The  inverse  problem  can  be 
simply  expressed  in  the  form  : 

St^  =  Gjt  •  (5m 

where 

X  indicates  one  of  {p,  g,  q,  a}. 

St^  is  an  AT -vector  for  corresponding  to  N  modes, 

G.  is  an  A/^  X  ^  Frechet  kernel  matrix  for  structural  inverse  prob¬ 
lem, 

<5m  is  a  model  correction  vector  for  k  unknowns. 

Keep  in  mind,  that  is  not  measurable  but  fortunately  GSDF  theory 
provides  a  way  to  compute  these  nonmeasurable  quantities  using  the 
mode  interference  relationship.  Thus,  in  actual  application  of  GSDF 
theory,  we  must  relate  to  ^t^  to  create  the  kernel  G^. 

From  (GJ.56-59),  we  can  transform  to  ^t^  as 

^tg  =  5rg{&f)  -1-  (1  -  -  Stgi&f)] 
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^tp  =  Stj^idif)  +  a[tcl  -dr^i&f)] 

Stf^  =  dTqf^Wf)-a8'Cg,{&f) 

where  a  =  (1  -  ^\){— — —).  Substituting  these  into  (GJ.84-87)  gives 

a>f 

the  following  equations  which  relate  the  four  GSDFs  to  inversion  ker¬ 
nels  Gx- 

Stp  =  1  •  C  •  la  H- 1  •  ^C(Gp  —  aGg)  —  S(Gq  —  aG^lJ  • 

Stq  =  1  •  S  •  la  -I- 1  •  ^C(Gq  —  aGa)  -I-  S(Gp  —  aGg)j  • 

St^  =  l-  (Cg  -  Sa)  •  la  +  1  •  S  •  1(1  - 

+  1  •  [-(Ca  +  Sg)(Gq  -  aGa)  +  (Cg  -  Sa)(Gp  -  aGg)]  •  Jm 

-i-l-[^fCGa  +  ^?SGg]-^m 

=  -  1  •  (Ca  +  Sg)  •  la  -I-  1  •  C  •  1(1  -  If  )^e 

+  1  •  [-(Cg  -  Sa)(Gq  -  aGa) "  (Ca  +  Sg)(Gp  -  aGg)]  •  <5m 

+  l-[-|fSGa  +  |fCGg]-^m 

Now,  we  have  a  general  form  of  inversion  problem  for  earth  structure 
in  terms  of  the  GSDF’s  5t^,  kernels  G^,  and  the  interference  effects 
C,S,Cx,Sx. 

At  this  point,  it  is  possible  to  outline  this  structure  inversion  prob¬ 
lem.  At  a  particular  frequency,  according  to  GSDF  theory,  using  cross¬ 
correlation  technique  we  can  obtain  four  measurements  as  our  "obser¬ 
vations"  (GJ.9-12)  and  evaluate  the  mode  interference,  which  we  will 
use  with  Gx  to  form  inversion  kernels.  Therefore,  we  perform  this 
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measiaring  procedure  at  several  frequencies,  we  then  have  the  "data 
vector"  and  "kernel  matrix"  ready  for  inversion.  In  the  following  sec¬ 
tions,  we  will  show  how  we  create  the  kernels  (  ). 


Kernel  Gp 


Kernel  Gp  relates  phase  velocity  changes  to  model  perturbations. 


^Tp  =  Gp  • 

The  n’th  component  in  the  differential  phase  velocity  vector  5x^,  SXp^ 
is  define  as 

R  R 


^^Pn  = 


^nobs  ^ 


nsyn 


R 


R  R 


Cn  + 


1-h(1-^  +  (^)2.-) 


R 


where 


R 

c„2  am 


•  (?m 


R  is  source-receiver  distance, 

Cn  is  phase  velocity  for  n’th  mode  of  synthetic  seismogram. 

In  the  matrix  form,  we  can  exactly  see  the  meaning  of  Stp^  as  the  total 
phase  delay  for  mode  1  when  model  m  changes  into  m  -i-  <5m. 


34 


r  R 

dci 

R 

dci 

R 

3ci 

dmi 

72 

dm2 

dms 

R 

dc2 

R 

dC2 

R 

dC2 

=  - 

r2 

dmi 

r2 

C2 

dm2 

r2 

C2 

dms 

Nxl 

.  . 

dcjq 

R 

dcjf 

R 

dCN 

R 

L  Off 

dmi 

7^ 

dm2 

r2 

dms 

(#1  mode) 
(#2  mode) 


5mi 

6  m2 

Sms 


(#N  mode) 


Smk 


\kxl 


^Nxk 


Kernel  Gg 


To  construct  kernel  Gg  which  relates  group  velocity  variation  to 
model  perturbations, 

STg  =  Gg  • 

the  n’th  component  in  differential  group  velocity  vector  Stg,  Szg^  is 
defined  as 


Sta  = 

gn  U 


R 


R 


nobs 


u 


nsyn 


R 


Un+SU, 


R 

Ul 


R 

uZ 


u„ 


t/. 


R 

U„ 


SUn  , 

Un  Uj 


R 


U 


su„  = 


2  n 


R  dU„ 


•  <5m 


.  am 

where 

Un  is  group  velocity  of  n’th  mode  synthetic  seismogram. 


In  matrix  form 
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''gN-^Nxl 


R 

3l7i 

R 

3Ui 

R 

dUi 

m 

dm-i 

U\ 

dm2 

U\ 

dm^ 

R 

3C72 

R 

dU2 

R 

dU2 

ui 

dm^ 

m 

dm2 

ui 

R 

dUr, 

R 

dUu 

R 

i)U„ 

dmi 

ul 

dm2 

ul 

dm^ 

(#1  mode) 

Smo 

(#2  mode)  _ 


i#N  mode) 


Kemel  Gq 

To  create  kernel  Gq  which  relates  attenuation  to  model  perturba¬ 
tions, 

(^Tq  =  Gq  • 

we  extend  (GJ.47)  into  a  general  matrix  form  and  reformulate  the 
above  equation  to  yield 


=  =?»[«-* -Q-iJ 


■q  2  f 


1  -  r  . 

2 


=  Gq 

Using  the  relationship  between  y  and  Q,  the  partial  derivative  of  y 
with  respect  to  model  perturbations  can  be  calculated. 

^Yn  ^  _  ^-1  (O 

3m  2co  3m  ”  2co  ^  3m 

Therefore,  the  partial  derivative  of  can  be  rewritten  as 
3m  3m  ”  2co  3m  co 
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The  w’th  row  (#  n  mode)  of  kernel  Gq  is 


(G,),  -  .,.(<»)  —  ^  +  2^7  to 


Gq  in  matrix  form  is 

.  fcpi  ^  ^gpi 

0)  dmi  2coj  dmi 

~  (co^  9/2  ,  Qi^  ^^02 

Tr.J - t: - r  - r - 


CO  drrii  dmi 


f  r^oi  I 

3^2  2coj  3ni2 

„  fC02  ^72  ^  Ql^ 

CO  dm2  2co2  ^”^2 


.  fcor,  ^Yn  ,  Qjv  ~  ((^Qn  ^Yn  ,  Qiv 

Tn.J  -  - r- -  -  T„„  -  :r - b  - - - 


0)  dmi  2cojy  dmi 


CO  dm2  2co^  dm2 


(#1  mode) 

(#2  mode) 

(,#N  mode) 

-(Nxk) 


Kernel  Ga 


The  last  kernel  is  Ga-  From  (GJ.48),  we  have  the  relationship 
between  Ga  and  Q. 

^  rg((o)[  -  Q  "M  =  Ga  • 

In  a  manner  similar  to  deriving  Gq,  we  can  use  the  Q  and  y  relation 
and  the  partial  derivatives  of  y  to  obtain  the  kernel  Ga  as 


(Ga)n 


^gn(^)  ^Qn 

2  dm 


Co„  ^  Q~n 

CO  dm  2  dm 


Ga  in  matrix  form  is 
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-  fcpi  M 

0)  dmi 

~  (^02  ^r2 

0)  dmi 


Ql^  ^Cq,  > 

2coi 

Qi^  9co2^ 

2co2  9/ni; 


cqi  ^ri  ^ 
CO  dm2 


V2 


^£o2 

CO  dm2 


+ 


Ql^  ^CQi^ 

2coi 

2co2  dm2^ 


(#1  mode) 
(#2  mode) 


~  ( ^Ojv  ^  Qjv 


Oiv 


CD  drui  2cq  drui 


( ^Yn 


CO  dm2  2co^ 


Q)J  3c"  ^ 


-Ojv 


9m 


2  y 


(#iV  mode) 


INxk) 


Partial  Derivatives  of  c,  t/,  and  y 

All  four  kernels  Gp  Gg  Gq  Ga  are  derived  in  terms  of  partial 
derivatives  of  phase  velocity  c,  group  velocity  U,  and  attenuation  y.  In 
this  section,  We  will  give  all  partial  derivatives  which  will  be  used  in 
creating  kernels. 

Perturbation  theory  is  used  to  obtain  phase  velocity  partials  with 
respect  to  medium  parameters,  and  a  numerical  method  which  intro¬ 
duced  by  Rodi  et  al.  (1975)  is  used  to  calculate  surface- wave  group- 
velocity  partial  derivatives.  The  following  equations  are  calculated  for 
a  single  mode  at  a  certain  frequency,  therefore  we  omit  the  subscript  n 
for  n’th  mode.  The  subscript  0  is  used  to  denote  the  value  before  intro¬ 
ducing  causal  Q.  The  subscript  v  can  be  substituted  by  P-wave  veloc¬ 
ity  a  or  S-wave  velocity  p.  pis  layer  density,  h  is  layer  thickness,  coj. 
is  a  reference  angular  frequency  used  for  introducing  causal  Q. 

Rayleigh  Wave 


1 ,  ,  O)  ^ 

c  =  Co  +  —  ln( — )  X 

JT  COy. 


^  5co  „ 


-1 

a 


r  = 


1  9co 


-1 

a 


38 


-•  .  Uo^,c-Co  _  2yUo 

u  =  Uq  1  +  (2 - )( - )  + - 

Cq  Cq  7UC0 


dc  dco^  1  6) 

dv  dv  nQv  cOr 


dc  dco  T dco  P  1  dcQ  a  1 

dp  dp  n  (Or  dp  2p  ^  da  2p 


dc  dco 


dc  1  (0  dco 

——  =  -  ln(— )  —  V 
dQ~^  n  (Or  dv 

dy  ^  (o  dcp  1  2/  dcp 

dv  2cl  dv  "  Co  dv 


dy  ®  Taco  ,  y?  dco  ,  a  1]  2y  dco 


iL.  =  __  J—'JO-I  +  1Z}L  (-  —)0-^  -  ZL.  Ziz 

dp  2cq  dp  2p  ^  da  2yO  “  Cq  dp 


dy  _  2y  dco 
dh  Cq  dh 

dy  _  (o  dcQ 
dQ-^  2cl  dv 


^^dU^  U  Up  ^c-Cq^  ^2yUQ  ^  dcp  q  ^  2  ^  1 

dv  dv  Uq  Cq  Co  TTCD  dv  Cq  Cq  Uq 


dv  Cq  Cq  dv  K(0 


^^^Uq  U  Uq  ^c-Cq^  ^2yUQ  ^  (Eo  )2  2  i- -  2 


dp  dp  \Uq  Cq  Cq 


dp  Cq  Cq  Uq 


+  ^11^(2-—)+  —  ^ 
dp  Cp  Cp  dp  7t(o 


U  Uo^c-Co^^2yUo 

dh  dh  Uq  Co  Cq  nco 

"bh  Co  Co  bh  nco 

w  tfp  t/p  9c  %ul  Sr 

3Qr'  <^o  «o  30;'  rria 


Love  Wave 


1  ,  .  CO  .  —  bCn  „„_i 
c  =  Co  +  -  ln( — )  X  3^  J^Q/ 
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U  =  Wo 
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dy  _  2y  dcQ 
dh  Cq  dk 

dy  ^  0)  dcp 

dQf  2cl  dfi  ^ 

^^dUo  U  Up  ^c-Cq^  ^2yUo  ^  dcp  q  ^  q  ^  1 

dfi  dp  Uq  Cq  Co  TTOJ  dp  Cq  Cq  Uq 

dp  Cq  Cq  dP  71(0 

^  =  ?£l  ^  ^0  (C-CQ)  ,  2/^0  ^  dCQ  ^Uq^2  0^0^  1 

dp  dp  Uq  Cq  Cq  TtCO  dp  Co  Co  Uq 

ic  U^,„_U^.dr^ 

dp  Cq  Cq  dp  71(0 

^^d^  U  Uq  ^c-Cq^  ^  2yUQ  ^  dcQ  ^Uq^2  q  ^  q  ^  1 

dh  dh  Uq  Cq  Cq  TTCO  dh  Cq  Cq  Uq 

dh  Cq  Cq  dh  TTCO 

dU  t/p  Up  dc  2ul  dy 

iQ't  Co  Co  dQ'^  ita  3Q}‘ 


3.4  Forming  Inversion  Kernel 

As  shown  above,  we  calculate  partial  derivatives  with  respect  to 
all  parameters  {a,  p,  p,  Qa,  Qp,  h)  instead  of  computing  partial  deriva¬ 
tives  for  shear  velocity  only.  The  intention  is  that  we  try  to  provide  all 
the  possible  tools  to  interpret  the  seismograms  as  completely  as  possi¬ 
ble,  and  as  automatically  as  possible.  Therefore,  it  is  the  users’ 
responsibility  to  choose  those  model  parameters  they  want  to  invert; 
and  only  the  chosen  part  will  be  used  to  assemble  the  inversion  kernel. 
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To  stabilize  the  inversion,  all  information  for  different  frequencies 
inside  the  inversion  kernel  are  weighted  according  to  their  frequency 
amplitudes.  This  weighting  procedure  greatly  improves  the  inversion 
stability. 

When  inverting  teleseismic  surface-waveforms,  sometimes  the 
body  waves  (e.g.  SS,  SSS)  or  some  unwanted  surface-waves  due  to 
improper  rotation  will  have  bad  influence  on  the  cross-correlation 
between  the  isolation  function  and  the  observed  seismogram,  and 
cause  signal  misalignment.  To  avoid  this  problem,  a  window  is  applied 
on  the  seismograms  before  doing  cross-correlation,  and  this  may  suc¬ 
cessfully  isolate  the  surface  wave  wavetrain  from  those  body  waves 
outside  the  surface  wave  wavetrain.  Although  prewindowing  proce¬ 
dure  is  applied,  the  same  trouble  still  may  happen  occasionally.  In 
such  a  situation,  the  information  at  that  frequency  must  be  rejected 
from  forming  inversion  kernel,  otherwise  it  will  plague  the  inversion 
for  its  strong  misaligned  "phase  delay"  or  "group  delay".  Figure  3.1 
shows  an  example  of  this  problem.  Due  to  improper  rotation,  the 
unwanted  Love-wave  waveform  appears  on  the  radial  component  seis¬ 
mogram  prior  to  the  Rayleigh  wave  wavetrain.  This  unwanted  Love 
wave  signal  causes  two  kind  mistakes,  wrong  identification  of  a  Gaus¬ 
sian  wavelet  (Figrue  3.1a)  and  signal  misalignment  (Figure  3.1b), 
which  can  not  be  incorporated  in  inversion. 

We  have  mentioned  that  users  have  the  power  to  decide  what 
parameters  will  be  inverted  for  during  the  inversion.  An  aspect  which 
strongly  related  to  this  decision  is  the  phrase  "using  what  kind  infor¬ 
mation  ?".  During  each  iteration  in  the  inversion,  we  have  calculated 
partial  derivatives  with  respect  to  parameters  for  each  layer  at  sev¬ 
eral  appointed  frequencies,  and  all  this  information  is  rearranged  to 
form  four  kind  delays  at  each  frequency.  After  manually  rejecting  mis¬ 
calculated  cross-correlations  at  some  frequencies,  users  have  to  decide 
which  combinations  of  Stp,  Stg,  St^,  Stg  they  will  use  to  invert  for  a 
particular  combination  of  model  parameters  (a,  P,  p,  Qa,  Qfi,  h)-  From 
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Fig.  3.1a.  This  is  an  example  showning  how  unexpected  signal  interference  affects 
the  extracting  procedure  in  GSDF  theory.  There  are  five  traces  presented  to  show  the 
different  processing  stages.  The  top  two  traces  are  the  prefiltered  isolation  filter  and 
observed  seismogram,  respectively,  The  third  trace  shows  the  Gaussian  filtered 
cross-correlation  at  a  target  period  of  20.0  seconds.  The  five  extracted  parameters  are 
shown.  The  dashed  curves  inside  the  envelope  are  from  the  synthetics  and  the  solid 
curve  are  from  the  observed  data.  The  fourth  trace  is  the  windowed  cross- 
correlagram  shown  in  the  third  trace.  The  bottom  trace  is  the  filtered  windowed 
cross-correlagram  and  its  five  Gaussian  wavelet  parameters  which  are  to  be  used  in 
further  processing.  Due  to  improper  rotation,  the  Love  wave  appears  on  radial  com¬ 
ponent  before  the  Rayleigh  wave  arrival.  The  Love  wave  wavetrain  may  causes  two 
kind  troubles  like  (a)  wrong  determination  of  Gaussian  wavelet  parameters,  and 
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easily  pull  the  model  close  enough  to  the  final  model  so  that  phase 
delays  can  be  used.  It  only  at  the  final  stage  when  the  synthetic  is  very 
close  to  the  observed  seismograms  that  the  phase  delay  information 
can  be  used  at  a  powerful  fine  tuning  of  the  model  since  there  is  no 
"cycle  skipping"  problem. 

3.5  Synthetic  Test 

A  simple  source  with  strike,  dip,  rake  angles  of  45°,  45°,  45°, 
respectively,  in  20  km  deep  was  used  in  this  synthetic  test.  Receiver  is 
located  at  1000  km  away  along  10°  in  azimuth.  The  "observed"  seismo¬ 
grams  were  created  for  a  two  layered  crust  with  an  upper  mantle  devi¬ 
ating  slightly  from  the  PREM  model  (Figure  3.2).  Here,  we  only  try  to 
invert  for  shear  velocity  structure,  so  all  the  other  parameters  related 
to  the  source  and  earth  model  are  assumed  known.  The  starting 
model  is  a  two-layered  crust  model  with  the  PREM  model  beneath  40 
km  depth.  It  is  clearly  seen,  in  Figure  3.3  that  the  starting  ssmthetic 
seismograms  are  very  far  away  "observed"  seismograms. 

After  12  iterations,  the  synthetic  seismograms  are  almost  the 
same  as  "observed"  seismograms  (Figure  3.4).  Checking  the  model  dif¬ 
ferences  between  the  "true"  model  and  the  final  model  (Figure  3.5),  we 
can  see  the  crust  structure  is  almost  matches  the  "true"  model,  except 
in  the  upper  mantle  where  the  surface  wave  doesn’t  provide  enough 
resolving  power  the  inverted  structure  causing  the  "zig-zag"  pattern. 
The  differences  in  Q  models  are  not  significant  at  these  frequencies 
and  this  distance. 

3.6  Real  Data  Test 

After  the  successful  synthetic  test,  we  wish  to  test  this  technique 
on  real  data.  The  April  14,  1995  Texas  earthquake  is  the  one  chosen. 
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Fig.  3.2.  The  starting  model  (dashed  line)  used  in  synthetic  test  of  GSDF  inversion 
algorithm  and  the 

According  to  Sipkin’s  USGS  solution,  the  epicenter  of  the  Texas  earth¬ 
quake  is  at  30.261°N  103.33°W  and  origin  time  is  00:32:55UT.  The 
current  collected  data  set  consists  of  89  broadband  stations  from  IRIS, 
Canadian  National  Seismological  Center  (CNSDC),  USGS,  UNAM, 
and  PASSCAL  experiment  (Figure  2.1). 

Some  source  parameters  were  reestimated  on  the  basis  of  fitting 
surface  wave  amplitude  spectra  (Table  2.2).  The  redetermined  source 
depth  is  23  km  deep,  with  strike,  dip,  rake  angles  of  114°,  64°,  -101°, 
respectively  with  M„,  =  5.6. 

3.7  Inversion  Procedure 

We  have  tried  three  different  inversion  runs,  the  first  two  exhibit 
some  significant  difficulties  in  matching  the  waveform  or  in  reason¬ 
ableness  of  the  resulting  model,  so  the  third  one  has  been  adapted  for 
inverting  structure.  In  performing  the  inversions,  we  use  an  earth 
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Fig.  3.3,  The  velocity  time  histories  of  'observed  seismograms'  (solid  line)  and  the 
those  from  the  starting  model  (dashed  line)  are  both  filtered  in  the  frequency  band 
0.01-0.05  Hz  by  using  a  Butterworth  filter  with  four  poles.  The  plotted  seismograms 
are  normalized  according  the  maximum  amplitude  of  each  component  in  current  fre¬ 
quency  band.  It  is  clear  to  see  that  the  starting  model  is  not  close  to  the  'true'  model 
and  this  may  demonstrate  the  ability  of  the  inversion  programs  to  resolve  structure. 
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Fig.  3.4.  After  12  iterations,  the  final  inversion  result  show  an  predicted  waveforms 
almost  identical  to  the  ’observed  seismograms’. 


flattening  approximation  to  use  plane-layered  surface-wave  theory  to 
generate  synthetics. 
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Fig.  3.5.  The  comparison  between  the  final  model  and  the  ’true  model’.  We  can  see 
that  the  2-layer  crust  is  very  close  to  the  ’true’  model,  but  that  in  the  upper  mantle 
the  model  shows  some  ’zig-zag’  pattern. 

The  first  run  consisted  of  a  joint  inversion  of  both  Rayleigh  and 
Love  wave  seismograms  for  shear  velocity  structure.  The  result  was 
that  that  the  synthetic  seismograms  tended  to  fit  the  largest  surface 
wave  amplitude,  and  ignored  the  small  amplitude  wavefields.  This 
arose  because  of  the  amplitude  level  weighting  used  to  stabilize  the 
inversion.  Therefore,  separate  inversions  for  Rayleigh  wave  and  Love 
wave  were  necessary. 

The  second  sequence  consisted  of  inverting  for  the  shear  velocity 
structure  from  the  Love  wave,  and  then  using  this  structure  as  an  a 
priori  shear  velocity  structure  so  the  the  Rayleigh  wave  provided  infor¬ 
mation  on  the  compressional  velocity  structure.  This  procedure  can 
match  both  Rayleigh  wave  and  Love  wave  waveforms  very  well,  but  we 
found  it  is  impossible  to  find  a  reasonable  explanation  for  the  anoma¬ 
lous  Poisson’s  ratios  in  the  inverted  model.  It  is  well  known  that  the 
surface  wave  are  insensitive  to  compressional  wave  velocity  structure, 
so  the  only  conclusion  for  this  is  that  the  shear  velocity  structure 
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inverted  from  Love  wave  is  not  adequate  for  Rayleigh  wave.  An  unac¬ 
counted  anisotropy  effect  may  cause  the  overcorrection  in  compression 
velocity  structure.  Figure  3.6  compares  waveforms  and  Figure  3.7 
shows  the  model  resulting  from  the  second  inversion  procedure. 
Although  the  synthetic  waveform  does  not  perfectly  fit  the  Rayleigh 
wave  (Figiue  3.6),  we  can  see  the  general  features  are  matched  using  a 
single  model  for  both  Love  and  Rayleigh  waves.  In  Figure  3.7  we  see 
that  the  P-wave  velocity  structure  has  lower  values  to  compensate  for 
the  high  shear  velocity,  and  vice  versa.  The  result  are  some  unexplain¬ 
able  Poisson’s  ratios  such  as  0.135  for  the  middle  crust  and  0.328  for 
the  lower  crust. 


p  P  a  Qb  Qa 


Fig.  3.7.  The  model  inverted  by  the  second  inversion  procedxire. 


Since  our  forward  S3mthetic  seismogram  algorithm  does  not 
include  anisotropy  effect,  we  modify  our  inversion  procedure  as  as  fol¬ 
lows:  we  invert  V p_sH  from  Love  wave,  and  another  separate  inversion 
of  V p_psv  from  Rayleigh  wave  on  the  vertical  component.  We  fixed  the 
Poisson’s  ratio  when  we  invert  V p  psv  from  Rayleigh  wave.  The  fixed 
Poisson’s  ratios  are  kept  the  same  as  the  original  input  model.  In  this 
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Fig.  3.6.  The  result  of  inversion  using  the  second  procedure.  This  inversion  procedure 
inverts  S  wave  velocity  from  Love  wave  and  after  that  inverts  P  wave  velocity  from 
Rayleigh  wave  by  assuming  no  anisotropy  effect.  This  shows  a  acceptable  waveform 
match  in  phase,  but  not  in  envelope. 


study,  we  set  the  Poisson’s  ratio  as  0. 25  for  crust,  0. 28  for  the  layers 
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between  40  km  and  220  km,  and  adopt  the  values  from  PREM  model 
for  those  layers  deeper  than  220  km. 

When  we  invert  for  shear-wave  velocity  structures,  the  attenua¬ 
tion  factors  can  be  determined  simultaneous  either  by  joint  inversion 
of  V p  and  Q  or  by  another  separated  Q  inversion.  The  Q  determina¬ 
tion  is  not  definitive,  since  we  have  a  lot  uncertainties  in  our  source 
and  velocity  structures  in  our  inversion,  so  the  only  objective  criteria 
for  determining  Q  is  the  envelope  shape  of  the  surface  wave  wavetrain. 
As  shown  on  Figure  3.8,  the  Rayleigh  wave  signal  at  large  distance  is 
well  dispersed  with  a  strong  Airy  phase.  This  Airy  phase  correspond¬ 
ing  is  affected  by  crustal  wave  propagation,  therefore,  the  envelope 
amplitude  of  the  Airy  phase  is  controlled  by  the  crustal  attenuation 
factors.  So  the  Q  structure  determined  is  sensitive  to  the  crust  and 
uppermost  mantle. 

3.8  Inversion  results 

Data  from  43  of  89  broadband  stations  were  inverted.  In  this 
report,  three  stations  were  selected  to  present  the  inversion  results 
and  the  waveform  fitting  success  will  be  shown  on  several  different 
frequency  bands. 

Station  ALE,  which  located  in  Arctic  region  is  the  farthest  station 
used  in  this  study.  There  are  some  interesting  features  for  the 
inverted  results.  Looking  at  the  waveform  fit  in  the  low  frequency 
range  (0.005-0.03  Hz  bandpass;  shown  as  Figure  3.8a),  it  is  clear  to  see 
that  the  synthetic  seismograms  successfully  match  the  observed  sur¬ 
face  wave.  Two  velocity  models,  PSV  and  SH  are  inverted  respectively 
from  Rayleigh  wave  and  Love  wave. 

For  ALE,  the  attenuation  factor  was  fixed  during  inversion.  A 
high  Q  structure  (Qg  =  1000)  is  used  for  crust,  a  low  Q  {Qs  =  100)  was 
adopted  for  structure  between  40  km  and  500  km,  and  using  Qs  =  143 
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Fig.  3.8a.  Waveform  fit  of  the  final  inverted  model  for  ALE  in  the  frequency  band  of 
0.005-0.03  Hz.  The  signal  which  arrives  at  1020  seconds  and  around  1300  seconds 
are  S  and  SS  phase,  respectively.  From  the  S  phase  waveform,  we  can  say  that  the 
source  time  function  used  in  this  inversion  is  a  little  short  but  is  close  enough. 


for  those  deeper  than  500  km.  And  from  Figures  3.8a,b,  the  synthetic 
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Fig.  3.8b.  (Cont’d).  Waveform  fit  in  the  0.005-0.05  Hz  frequency  band. 


Rayleigh  wave  and  Love  wave  amplitudes  only  show  small  deviation 
from  observed  seismogram,  therefore  the  Q  model  is  considered  ade¬ 
quate.  We  also  notice  this  sharp  Q  contrast  between  crust  and  mantle 
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Fig.  3.9.  The  inverted  models  for  ALE.  Two  models  (PSV  and  SH)  were  inverted  for 
Love  wave  and  Rayleigh  Wave.  From  these  two  models,  there  is  one  slight 
anisotropic  zone  above  200  km. 

is  a  common  character  for  those  stations  located  inside  the  North 
America  craton  region. 

The  second  station  is  FCC  which  located  on  the  west  shoreline  of 
Hudson  Bay  and  is  in  the  center  of  the  North  America  craton.  The 
inverted  models  (Figure  3.10)  show  high  shear  velocity  for  upper  man¬ 
tle  but  not  as  high  as  SNA  model  (Grand  and  Helmberger,  1984).  The 
current  inversion  result  shows  4.5%  anisotropy  between  70  km  and 
140  km.  The  synthetics  fit  data  well  at  low  frequency  (Figure  3.11a) 
and  can  fit  the  fundamental  mode  as  high  as  0.1  Hz  (Figure  3.11),  but 
have  difficulty  to  produce  some  higher  mode  arrivals  at  800  and  940. 
the  time  range. 

The  final  station  is  PAS.  The  wave  propagates  through  the  south¬ 
ern  Rocky  Mountains.  The  inverted  model  (Figure  3.12)  does  not 
require  anisotropy  and  the  shear  velocity  model  is  very  close  to  the 
TNA  model  (Grand  and  Helmberger,  1984).  The  synthetics  fit  the  S 
phase  which  arrives  around  350  seconds  and  surface  wave  very  well 
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Fig.  3.10.  The  final  inverted  model  for  FCC.  There  is  4.5%  anisotropy  effect  exists 
between  70  and  140  km. 

(Figures  3.13abc)  and  this  suggests  that  there  is  a  velocity  discontinu¬ 
ity  at  220  km  which  can  not  be  seen  in  the  TNA  model.  The  Q  is  very 
low,  with  the  average  Q  for  crust  of  lower  than  200.  However  we  found 
one  interesting  feature  about  the  Q  behavior.  The  Q  values  between 
40  and  220  km  are  slightly  higher  than  PREM  model,  i.e.,  a  low  Q 
crust  underlain  a  slightly  high  Q  upper  mantle  with  respect  to  the  ref¬ 
erence  model. 

3.9  Conclusion 

We  implemented  the  (generalized  Seismological  Data  Functionals 
technique  of  Gee  and  Jordan  (1992)  in  a  surface-wave  waveform 
inversion  algorithm.  A  simple  synthetic  test  shows  its  robust  inver¬ 
sion  ability.  After  a  successful  S5mthetic  test,  we  useed  this  inversion 
algorithm  on  real  data  to  further  test  its  ability.  The  Texas  earth¬ 
quake  (30.26  °N  103.33°W,  00:32:55UT,  April  14,  1995)  is  a  very  good 
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AZ=9.879  BAZ=196.549  DIST=  3244.330 

Fig.  3.11a.  The  waveform  fitting  for  final  model  in  FCC  are  shown  at  three  frequency 
bands  :  (a)  0.01-0.03  (b)  0.01-0.05  (c)  0.01-0.1  Hz.  The  SS  phase  arrives  at  660  sec¬ 
onds.  The  inverted  model  can  fit  fundamental  mode  Love  wave  and  Rayleigh  wave 
waveforms  as  high  as  0.1  Hz,  but  it  lacks  the  ability  to  simulate  the  higher  modes. 


earthquake  for  this  purpose  because  it  was  well  recorded,  its  source 
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Fig.  3.11b.  (Cont’d).  (b)  0.01-0.05  Hz. 


depth  is  well  constrained,  the  focal  mechanism  is  fairly  determined, 
and  the  seismic  moment  is  constrained  by  long  period  surface  waves. 
The  inversion  results  are  excellent  and  show  some  interesting 
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Fig.  3.11c.  (Cont’d).  (c)  0.01-0.1  Hz.  The  SS  phase  arrives  at  660  seconds.  The 
inverted  model  can  fit  fundamental  mode  Love  wave  and  Rayleigh  wave  waveforms 
as  high  as  0.1  Hz,  but  it  lacks  the  ability  to  simulate  the  higher  modes. 


features.  For  the  craton  there  is  some  evidence  for  anisotropy  and 
crustal  Q  is  high.  For  the  mountain  region,  although  the  inverted 
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Fig.  3.12.  The  inverted  models  for  PAS.  There  is  no  clear  anisotropy  effects.  The 
crustal  Q  is  very  low. 

model  show  very  similar  shear  velocity  structure  like  TNA  model,  but 
from  the  S-wave  waveform  the  model  prefers  a  velocity  discontinuity 
at  220  km.  These  features  worth  more  effect  on  them. 
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Fig.  3.13b.  (Cont’d).  (b)  0.01-0.05  Hz. 
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Fig.  3.13c.  (Cont’d).  (c)  0.01-0.1  Hz. 


CHAPTER  4 

GENETIC  ALGORITHMS 


Most  of  seismological  inverse  problems  are  nonlinear.  The  tech¬ 
niques  used  to  solve  such  nonlinear  problems  can  be  concluded  as  two 
groups.  The  strategy  of  first  group  is  to  hnearize  nonlinear  problems, 
then  use  iterative  processes  to  seek  a  better  solution  by  using  the  gra¬ 
dient  information  of  misfit  function.  The  strategy  of  second  group  is  to 
directly  search  the  model  space,  and  find  the  acceptable  models. 

The  methods  such  as  least  squares  method,  steepest  descent,  and 
conjugate  gradients  belong  to  the  first  group.  Although  these  methods 
are  widely  used  in  seismology,  there  is  a  well-known  disadvantage  that 
is  the  requirement  of  a  good  starting  model.  For  studying  large  scale 
low-frequency  structures  or  deep  features  which  beneath  the  litho¬ 
sphere,  this  would  not  be  a  real  problem  because  the  researches  of  the 
past  half  century  already  provide  some  good  starting  models  such  as 
PREM,  IASPEI91  (Kennett  and  Engdahl,  1991).  In  the  other  hand,  for 
stud3dng  the  crustal  or  lithospheric  structures,  a  good  starting  model 
may  not  be  available  because  the  crust  or  lithosphere  is  the  most 
structural  heterogeneous  region  in  the  Earth.  However,  to  study  the 
lithosphere  is  very  important  because  the  evolution  history  of  litho¬ 
sphere  is  hidden  in  its  structure.  We  have  to  investigate  the  structure 
of  the  crust  and  lithosphere  all  over  the  world  to  know  the  evolution 
history,  to  understand  the  tectonic  processes,  to  make  better  prediction 
of  seismic  activities,  and  to  lower  the  damage  of  hazards. 

Base  on  such  a  situation,  we  need  to  find  a  way  to  investigate 
many  possible  models,  i.e.  to  search  the  whole  model  space  and  find 
out  the  possible  models  of  structure  and  their  variations.  The  methods 
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like  Monte  Carlo  method,  simulated  annealing  (SA),  and  genetic  algo¬ 
rithms  (GA)  are  belong  to  this  group.  The  Monte  Carlo  method  is  a 
random  search  method  which  has  been  used  in  seismology  for  a  long 
time  (e.g.  Press,  1968;  Keilis-Borok  and  Yanovskaja,  1967).  Press 
(1968)  showed  a  successful  experiment  which  use  Monte  Carlo  method 
to  search  the  model  that  can  produce  correct  body-wave  travel-times, 
sxirface-wave  dispersion,  the  earth’s  free  oscillation  periods,  mass,  and 
moment  of  inertia.  Six  models  were  found  from  about  five  million  ran¬ 
domly  generated  models.  As  pointed  out  by  Press  (1968),  the  reason 
that  using  Monte  Carlo  method  is  because  it  offers  the  advantage  of 
exploring  the  range  of  possible  solutions  and  indicate  the  degree  of 
uniqueness  achievable  with  currently  available  geophysical  data. 
Examine  these  six  models,  we  can  find  that  the  structure  for  the  lower 
mantle  are  pretty  consistent,  but  large  variation  in  the  upper  mantle. 
However,  Monte  Carlo  method  has  its  own  disadvantage.  As  pointed 
out  by  Keilis-Borok  and  Yanovskaja  (1967),  Monte  Carlo  method  does 
not  use  information  obtaining  from  previous  trials  in  the  next  trial. 

For  recent  years,  simulated  annealing  and  genetic  algorithms 
methods  are  very  popular  in  seismological  society.  The  simulated 
annealing  method  is  stimulated  from  the  crystalizing  process  observed 
in  chemistry.  The  genetic  algorithm  is  inspired  from  the  evolution  pro¬ 
cess  observed  in  biological  science.  These  two  algorithms  are  better 
search  methods  than  Monte  Carlo  method  because  they  use  the  infor¬ 
mation  obtained  in  previous  trials.  The  seismological  applications 
such  as  estimation  of  residual  statics  (Rothman,  1985  and  1986),  wave¬ 
form  inversion  of  reflection  data  (Sen  and  Stoffa,  1992;  Stoffa  and  Sen, 
1991;  Sambridge  and  Drijkoningen,  1992),  earthquake  h3rpocenter 
location  determination  (Sambridge  and  Gallagher,  1993),  and  receiver 
function  inversion  (Shibutani,  Sambridge,  and  Kennett,  1996).  There 
is  one  comon  thing  for  all  these  applications,  that  is  they  all  deal  with 
nonlinear  problems,  and  particularly  in  some  cases  which  due  to  com¬ 
plexity  of  structure. 
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For  example,  in  the  receiver  function  inversion,  Ammon  et  al. 
(1990)  showed  that  the  final  models  were  dependent  on  the  initial 
models.  Shibutani  et  al.  (1996)  showed  that  using  genetic  algorithm 
can  estimate  an  average  model  which  is  more  stable  and  less  depen¬ 
dent  on  the  starting  assumptions. 

We  can  see  that  both  simulated  annealing  and  genetic  algorithms 
are  good  ways  to  perform  the  imcertainty  assessment  in  a  complicated 
nonlinear  problem.  As  stated  by  Sambridge  and  Drijkoningen  (1992) 
on  SA  and  GA  methods:  "any  problem  feasible  by  one  could  also  be 
tackled  by  the  other".  To  select  which  methods  (GA  or  SA)  is  better  on 
surface-waveform  modeling  problem,  we  need  to  know  both  algorithms 
and  the  purpose  of  our  application.  As  is  well  known,  generating 
multi-mode  surface  wave  synthetics  is  computationally  intensive.  So 
the  computation  time  will  be  a  crucial  factor  for  selection. 

4.1  Workflow  of  Simulated  Annealing  method 

The  computation  procedures  of  simulated  annealing  is  like  follow¬ 
ing. 

•  start  from  an  arbitrary  model 

•  temperature-loop  :  at  temperature  T  -T^-k  -  ST 

□  parameter-loop  :  for  model  parameter  <Sj,  i  =  1,  •  •  •,  m 

fix  all  other  parameter  value  except  Si 

for  the  parameter  Si,  there  is  n  possible  values. 

o  possible-value-loop  :  <Sy ,  j  =  l,  -,n 

calculate  the  energy  function  EiSy)  which  is  the  nor¬ 
malized  cross-correlation  of  observed  and  synthetic 
seismogram. 

o  end  of  possible-value-loop 

□  end  of  parameter-loop 
compute  the  probability  distribution 
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P(Sy)  = 


r-EiSij) 

exp[  ^  ] 

^  —E(Sii) 

X  exp[^r^] 
;=i  ^ 


•  end  of  temperature-loop 

We  can  see  that  for  each  temperature,  it  is  necessary  to  perform  {m  ■  n) 
forward  computation  of  S5mthetics.  If  there  is  k  step  in  lowering  tem¬ 
perature  to  reach  the  global  minimum,  the  total  forward  computation 
will  be  (^  •  m  •  n).  The  problem  is  that  there  is  no  rule  in  choosing  the 
starting  temperature  Tq  and  increment  of  temperature  difference  5T. 
Basu  and  Frazer  (1990)  designed  a  sequence  test  runs  to  find  the  criti¬ 
cal  temperatvu-e.  Even  though,  it  is  still  too  time  consuming  for  sur¬ 
face-waveform  modeling. 


4.2  Workflow  of  genetic  algorithm 


For  a  m  member  society  evolving  through  n  generation,  the  com¬ 
putation  sequence  of  genetic  algorithm  is  as  following. 

•  Randomly  generate  m  individuals  as  first  generation. 


Generation-loop  :  for  generation  =  1,  •  •  •,  n 

Compute  m  synthetic  seismograms  (individuals) 

Evaluate  each  individual’s  performance;  i.e.  calculate  the 
goodness-of-fit. 

□  population-loop  :  for  child  =  1,  •  •  •,  m 

Base  on  the  individuals’  performance  (probability) 
to  select  them  as  parents; 
to  change  parents’  DNA;  and 
to  apply  possible  mutation  on  their  children. 

□  end  of  population-loop 


end  of  generation-loop 
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We  can  see  it  is  possible  to  perform  GA  on  a  small  population  society. 
This  will  be  more  computational  efficient  than  SA  in  the  surface- 
waveform  modeling  problem.  However,  perform  GA  on  a  small  popula¬ 
tion  society  has  its  own  risk,  as  pointed  out  by  Sambridge  and  Dri- 
jkoningen  (1992).  That  is  when  society’s  members  are  not  close  to 
global  solution,  the  relatively  good  individual  in  the  society  will  multi¬ 
ple  itself  and  dominate  the  population  (i.e.  trapped  in  local  minimum). 
This  problem  can  be  solved  by  increase  the  population  size.  However, 
our  purpose  for  using  GA  in  surface-waveform  modeling  is  not  relying 
on  GA  to  reach  the  global  minimum.  Instead,  we  prefer  to  have  sev¬ 
eral  runs  to  see  the  possible  uncertainty,  get  the  rough  idea  about  the 
structure,  and  find  some  good  initial  models  for  other  inversion  tech¬ 
niques.  In  our  test,  a  GA  run  will  take  1  to  3  hours  of  CPU  time  in 
SUN  ULTRA  1  (167MHz)  workstation.  It  is  affordable  to  have  several 
reruns  if  we  find  it  trapped  in  such  situation.  So  we  choose  the  genetic 
algorithm  for  modeling  surface-waveform.  In  the  following  sections, 
we  will  address  several  technical  issues  in  appling  GA  in  surface- 
waveform  modeling. 

4.3  Smoothing  mechanism 

Since  GA  is  one  type  of  random  search  method,  there  is  no  strong 
constraint  between  parameters  (layer  S-velocity  in  this  case).  Usually, 
there  will  be  a  strong  zig-zag  pattern  in  the  velocity  model;  this  is 
undesirable  for  the  purpose  in  seismology  if  the  emphasis  is  on  deter¬ 
mining  the  simplest  acceptable  model.  To  reduce  this  pattern  in  veloc¬ 
ity  model,  we  smooth  the  layer  velocity  by  considering  the  adjacent 
velocity  contrast  without  change  the  vertical  travel  times.  As  shown 
on  Figure  4.1,  the  original  model  (solid  line)  has  a  strong  zig-zag  pat¬ 
tern  but  the  smoothing  mechanism  reduces  the  contrast  between  layer 
shear  velocities  (dashed  line  represents  the  after  smoothed  model).  We 
are  thus  attempting  to  find  the  smoothest  model  consistent  with  data. 
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Fig.  4.1.  In  GA,  models  are  randomly  generated,  so  there  always  have  some  ’ZIG¬ 
ZAG’  patterns.  For  obtaining  a  smooth  background  velocity  model,  a  smoothing 
mechanism  is  introduced.  Without  changing  the  vertical  shear  wave  travel  time,  the 
smoothed  model  (dashed  line)  has  less  ’ZIG-ZAG’  pattern  than  the  original  model 
(solid  line). 


The  introducing  of  this  smoothing  mechanism  can  be  viewed  as  a 
heavy  mutation  case  in  GA.  Of  course,  the  best  model  may  not  survive 
under  such  mutation  through  final.  But,  we  can  find  it  from  the  record 
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of  models  of  each  generation.  This  may  help  to  escape  from  some  local 
TTiim'Trmm  in  some  cases  and  provide  driving  force  for  evolution. 

4.4  Generation  number  and  population  size 

Although  GA  is  a  global  search  method  which  can  potentially  find 
the  global  minimum,  we  did  not  set  that  as  our  goal  in  this  study.  Due 
to  the  intensive  computation  load  of  generating  multi-mode  surface- 
wave  synthetics,  we  limited  our  computations  to  a  small  population 
size  and  only  perform  it  through  finite  generations.  We  hoped,  via 
using  GA  search  method,  to  get  some  good  models  to  use  as  starting 
models  for  other  inversion  algorithms.  To  understand  what  generation 
number  is  sufficient  for  our  purpose,  we  have  tested  the  consequences 
of  a  large  generation  number.  In  this  test,  shown  as  in  Figure  4.2,  we 
performed  500  generations  and  find  that  after  50  generations  model 
improvement  proceeds  less  rapidly,  indicating  a  degree  of  convergence. 
Therefore,  in  the  subsequent  tests,  we  will  only  perform  50  genera¬ 
tions  for  a  small  population  (i.e.  20),  and  this  will  only  consume  1  to  3 
hours  of  CPU  time  in  SUN  ULTRA  1  workstation. 

In  some  cases,  we  do  find  the  GA  trapped  in  a  local  minimum. 
Usually  this  situation  is  associated  with  other  difficulties  like  cycle¬ 
skipping  problem  for  teleseismic  waveforms.  This  will  be  discussed 
later. 

4.5  Criteria  of  goodness-of-fit 

A  surface-wave  signal  has  a  longer  duration  and  a  more  compli¬ 
cated  waveform  behavior  than  any  single,  pulse-like  body-wave  phase. 
To  model  such  long-duration  complicated  waveforms,  there  is  a  cycle¬ 
skipping  problem  which  may  produce  an  unreasonably  low  or  high 
velocity  model.  In  addition  when  processing  surface-wave  data,  we 
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Fig.  4.2.  To  know  how  many  generations  is  necessary  for  surface-waveform  modeling, 
an  experiment  with  large  generation  number  (500)  is  tested.  The  result  shows  that 
GA  can  find  a  fairly  good  model  within  50  generations.  After  that,  the  model 
improvement  proceeds  less  rapidly. 

cannot  shift  the  synthetic  seismogram  to  match  the  observed  arrival 
time,  a  well  adapted  technique  in  processing  body  wave  data  such  as 
receiver  function  inversion.  Due  to  this  effect,  an  L2  norm,  such  as 
used  by  Gomberg  and  Masters  (1988),  may  not  be  suitable  for  quanti¬ 
fying  surface-waveform  goodness-of-fit.  Instead  we  choose  a  cross¬ 
correlation  as  our  criteria  of  goodness-of-fit  to  circumvent  the  oscilla¬ 
tory  signal  character,  and  to  focus  on  agreement  of  waveform  shapes. 

Surface-waves  usually  have  a  broad  frequency  content,  which 
means  a  single  cross-correlation  measurement  only  represents  the  fit 
of  the  largest  amplitudes,  which  are  typically  high-frequency  for 
crustal  earthquakes.  This  will  only  resolve  the  very  shallow  part  of 
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structure  but  leave  the  deeper  structure  uncertain  with  high  variation. 
To  overcome  this  problem,  we  divide  the  frequency  range  of  interest 
into  several  subranges  and  evaluate  the  cross-correlation  of  narrow 
band-filtered  observed  and  S3mthetic  seismograms  for  each  subrange. 
An  averaged  cross-correlation  value  of  these  subrange  cross¬ 
correlations  is  used  as  our  goodness-of-fit.  For  example,  for  ANMO,  we 
divided  the  period  range  (10-50  sec)  into  4  intervals  ;  (10-20  sec), 
(20-30  sec),  (30-40  sec),  (40-50  sec).  Using  these  period  intervals  as  the 
ranges  for  bandpass  filtering  observed  and  S3nithetic  seismograms,  a 
cross-correlation  value  is  computed  for  each  interval  and  an  average 
cross-correlation  is  used  as  our  goodness-of-fit. 

4.6  Test  on  the  western  Texas  earthquake 

We  apply  this  GA  search  method  to  the  April  14,  1995  Texas 
earthquake  (30.26  °N  103.33°W,  00:32:55UT).  The  source  depth  of  the 
Texas  event  is  23  km,  with  strike,  dip,  rake  angles  of  114°,  64°,  -101°, 
respectively  with  Mj^,  =  5.6. 

Three  stations  (ANMO,  TUG,  WMOK)  are  selected  to  show  the 
ability  of  the  GA  search  method.  For  each  station,  three  plots  of  final 
good  models,  waveform  fitting,  and  cross-correlation  of  different  period 
ranges  of  the  best  result  are  shown.  The  first  plot  shows  the  search 
bounds  (thickest  dashed  line),  the  good  models  (thin  dashed  line) 
which  have  goodness-of-fit  greater  than  a  certain  value,  and  the  best 
model  (black  solid  line).  The  second  plot  shows  the  waveform  fit  of  the 
best  model.  The  observed  seismogram  is  drawn  as  a  black  line  and  the 
synthetics  as  a  thin  dashed  line.  The  third  plot  shows  the  cross¬ 
correlation  measurements  of  the  best  searched  model  in  the  different 
period  ranges  which  gives  us  the  idea  how  good  can  this  model  fits  the 
data.  We  also  can  see  the  reason  why  using  averaged  cross-correlation 
as  our  criteria  of  goodness-of-fit  because  the  broadband  (e.g.  10-50  sec) 
waveform  at  top  row  is  dominated  by  high  frequency  signal  such  as 
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waveform  on  the  second  row  (e.g.  10-20  sec). 
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Fig.  4.3b.  (cont’d).  (b)  Searched  models  for  Rayleigh  wave. 

Examining  the  results  for  ANMO  (Figures  4.3,  4.4,  4.5),  we  found 
that  our  best  final  model  for  Rayleigh  wave  has  a  unreasonably  low 
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Fig.  4.3a.  Using  GA  to  model  the  surface-waveform  of  station  ANMO.  In  this  test, 
two  GA  runs  were  conducted  for  Love  and  Rayleigh  wave  respectively.  The  searched 
models  are  shown  in  these  figures,  the  heavy  lines  are  the  search  bounds,  the  thin 
lines  are  searched  model  which  have  their  goodness-of-fit  greater  than  a  certain 
value  shown  at  the  left  bottom,  and  the  best  model  is  plotted  as  thick  black  line,  (a) 
Searched  models  for  Love  wave. 
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FILTER  Flo=0.010  Fhi=0.100  (Hz)  Norder=  4 

LON=-103.327  LAT=30.261  DEPTH=  23.00 

AZ=331.340  BAZ=149.660  DIST=  596.790 


Fig.  4.4.  For  ANMO,  the  observed  waveforms  (solid  line)  and  predicted  waveforms 
(dashed  line)  which  generated  for  the  best  searched  models  for  Love  and  Rayleigh 
wave  at  frequency  range  0.01  to  0.1  Hz. 


velocities  for  layers  deeper  than  50  km.  Also  on  the  cross-correlation 
diagram,  we  can  see  that  at  the  40-50  second  period,  the  envelope’s 
maximum  of  the  Z  component  is  off  central  position.  This  is  the  flaw  of 
currently  used  criteria  which  can  not  overcome  the  cycle-skipping 
problem.  Combining  the  correlation  coefficient  at  zero  lag  with  the  lag 
shift  of  the  maximum  correlation  may  be  another  goodness  of  fit 
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Fig.  4.5.  For  ANMO,  the  filtered  cross-correlations  of  observed  and  synthetic  wave¬ 
forms  at  different  frequency  bands.  The  number  at  the  right  of  cross-correlation 
traces  is  the  cross-correlation  value  at  zero-lag  which  is  used  to  construct  the  value 
of  goodness-of-fit. 

criterion  to  use  in  the  future.  However,  compare  two  best  models 
obtaining  from  separated  GA  searchs  for  Love  and  Rayleigh  wave,  we 
can  find  these  two  models  show  a  very  similar  model  for  upper  50  km. 
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Fig.  4.6a.  Separated  GA  searchs  for  Love  and  Rayleigh  wave  recorded  at  TUC.  (a) 
Searched  models  for  Love  wave. 


At  TUC  (Figures  4.6,  4.7,  4.8,  and  4.9),  we  tested  this  search 
scheme  on  three  cases:  using  Rayleigh  wave  only,  using  Love  wave 
only,  and  using  both  Rayleigh  and  Love  wave.  The  final  results  show 
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GOODNESS  >  0.8  00  BEST=  0.930 

Fig.  4.6b.  (cont’d).  (b)  Searched  models  for  Rayleigh  wave. 

that  waveform  fitting  from  the  separated  search  are  better  than  those 
from  the  joint  search.  But  in  the  macroscopic  view,  they  all  have  a 
very  similar  velocity  gradient  in  the  crust.  This  may  illustrate  that 
the  crustal  structure  is  well  resolved.  In  the  uppermost  mantle, 
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Fig.  4.7.  A  joint  GA  search  using  both  Love  and  Rayleigh  waveforms  at  TUC. 

Figure  4.6(a)  shows  that  a  fixed  half-space  beneath  80  km  may  not  be 
appropriate  and  cause  the  model  which  between  60  to  80  km  to  lower 
their  values  to  compensate  this  high  velocity  half-space.  However, 
from  Figure  4.6(b),  the  Rayleigh  wave  seems  prefer  this  high  velocity 
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Fig.  4.8a.  For  TUC,  filtered  waveforms  at  fi*equency  range  0.01  to  0.1  Hz  for  observed 
and  predicted  seismograms,  (a)  The  predicted  seismograms  were  generated  using  the 
best  models  from  sepeirated  searched  for  Love  and  Rayleigh  wave  (Figures  4.6ab). 


half-space.  So  we  need  to  conduct  more  tests  to  see  if  there  exists  an 
anisotropy  zone  beneath  the  propagation  path. 

At  WMOK  (Figures  4.10,  4.11,  and  4.12),  the  waveform  of  the  best 
model  fits  the  observed  seismogram  very  well,  not  only  for  the  funda¬ 
mental  mode  but  also  for  the  first  higher  mode  Rayleigh  wave.  The 
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Fig.  4.8b.  (cont’d).  (b)  The  predicted  seismograms  were  generated  using  the  best 
model  from  a  joint  searched  for  both  Love  and  Rayleigh  wave  (Figure  4.7). 


propagation  path  through  the  west  Texas  region  is  only  655  km.  This 
region  under  the  propagation  path  is  a  uniform  platform  between 
Rocky  mountain  and  Ouachita  orogenic  belt.  From  the  final  model,  we 
can  see  the  existence  of  a  transitional  crust-mantle-boundary  between 
35  and  50  km. 
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Fig.  4.9a.  For  TUC,  the  cross-correlations  for  the  observed  and  S3nathetic  seismo¬ 
grams  which  were  generated  using  GA  searched  models,  (a)  The  synthetics  were 
computed  using  the  best  models  from  separated  searched  models  for  Love  and 
Rayleigh  wave  (Figures  4.6ab). 
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Fig.  4.9b.  (cont'd).  (b)  The  synthetics  were  computed  using  the  best  model  from  a 
joint  search  for  both  Love  and  Rayleigh  wave  (Figure  4.7). 


4.7  Test  for  teleseismic  traces 


In  our  test,  we  found  that  our  implement  of  GA  described  above  is 
only  work  for  short  epicenter  distance  records.  For  the  teleseismic 
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Fig.  4.10.  The  models  from  GA  search  for  Rayleigh  wave  recorded  at  WMOK. 

records,  this  algorithm  doesn’t  work  simply  because  the  cycle-skipping 
problem.  However  the  surface  wave  at  teleseismic  distance  samples 
the  upper  mantle,  and  is  very  important  to  provide  information  and 
constraints  on  the  upper  mantle.  We  need  to  find  a  way  to  apply  GA 
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Fig.  4.11.  The  waveform  comparison  of  observed  and  synthetic  seismograms  for 
WMOK  (Figure  4.10). 

search  technique  to  teleseismic  records. 

To  apply  GA  search  technique  to  teleseismic  seismograms,  we 
change  the  way  of  parameterization.  We  use  three  20  km  thick  layers 
over  a  half-space,  which  corresponding  to  the  upper  crust,  lower  crust, 
uppermost  mantle,  and  the  upper  mantle.  From  the  reports  on  the 
structure  of  lithosphere,  we  notice  that  the  average  velocity  for  the 
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Fig.  4.12.  The  cross-correlation  of  the  best  GA  searched  model  for  WMOK  (Figure 
4.10). 

uppermost  mantle  (40-60  km)  is  greater  than  the  lower  crust  (20-40 
km),  and  the  lower  crust  velocity  is  higher  than  upper  crust  (0-20  km). 
In  the  upper  mantle  (half-space  in  this  parameterization),  the  average 
velocity  should  not  deviate  too  much  from  the  uppermost  mantle.  So 
we  place  these  physical  constraints  on  our  GA  searched  models,  i.e.  the 
velocities  of  the  first  three  layer  should  be  monotonic  increasing  and 
the  difference  between  the  uppermost  mantle  and  upper  mantle  should 
not  greater  than  0.2  km/sec. 

The  test  results  show  a  good  improvement  in  searched  model,  it 
overcomes  the  cycle-skipping  problem.  Three  results  (HRV,  FRB, 
LMN)  are  shown  here  to  show  the  successful  GA  searched  results. 
However,  this  strategy  has  its  own  weakness.  We  will  use  the  GA 
searched  results  for  INK  to  illustrate  the  problem. 

In  HRV  (Figures  4.13,  4.14),  the  best  GA  searched  model  can  pre¬ 
dict  the  waveform  very  well  at  the  frequency  range  0.01  to  0.05  Hz. 
The  models  show  a  fairly  small  uncertainty  and  the  best  model  is 
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consistent  with  the  velocity  structure  of  stable  continents.  In  FRB 
(Figures  4.15,  4.16),  the  GA  search  show  its  ability  to  find  the  models 
to  predict  the  well  dispersed  waveform.  We  also  notice  the  envelope  of 
synthetics  is  small  than  observed  seismograms  which  may  reflect  that 
the  model  is  not  good  enough  to  produce  correct  amplitude.  From  Fig¬ 
ure  4.16,  we  can  see  that  lower  crust  velocity  may  be  too  low  and 
uppermost  mantle  velocity  may  be  a  little  high.  In  LMN  (Figures  4.17, 
4,18),  the  waveform  fit  is  all  right  and  its  model  looks  reasonable,  but 
we  can  see  the  synthetics  does  not  match  the  observed  Airy  phase. 
This  may  be  caused  by  the  criteria  of  goodness-of-fit.  In  appl3dng  GA 
search  technique  to  teleseismic  seismograms,  we  divided  the  period 
range  10-70  seconds  into  6  subranges.  The  period  content  of  wave 
traveling  through  crust  is  mainly  shorter  than  30  seconds,  therefore 
the  criteria  of  goodness-of-fit  has  more  weighting  for  low-frequency  sig¬ 
nals  and  may  not  properly  represent  the  crust  wave.  The  problem  is 
more  clear  in  INK  (Figures  4.19,  4.20).  Wen  can  see  that  due  to  the 
lack  of  proper  representation  of  crustal  signals,  GA  search  technique 
fail  in  searching  good  model  for  the  noisy  traces  like  INK  which  suffer 
from  multipathing  problem. 

4.8  DISCUSSION 

From  this  experiment,  we  found  that  GA  search  method  works 
well  for  regional  seismograms  (A  <  1000  km)  but  not  good  for  teleseis¬ 
mic  traces  (A  >  2500  km  ).  The  reasons  for  this  are  as  follow; 

•  First,  the  regional  seismograms  have  a  more  concentrated  energy 
envelope  instead  of  a  well-dispersed  wavetrain  seen  in  teleseismic 
surface  waves.  Therefore,  the  cycle-skipping  problem  is  less 
severe  in  this  case. 

•  Second,  the  goodness-of-fit  criteria  may  be  too  simple.  The  crite¬ 
ria  only  utilize  the  amplitude  information  of  the  cross-correlation 
but  ignore  the  time-shift  information. 
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Fig.  4.13.  The  waveform  fit  of  the  best  GA  searched  model  for  HRV. 


A  different  parameterization  which  designed  for  teleseismic  seis¬ 
mograms  is  working  well  for  some  high  S/N  ratio  with  strong  Airy 
phase  traces. 

Another  important  result  of  genetic  algorithms  is  that  they  pro¬ 
vide  a  suite  of  possible  models.  The  distribution  of  possible  solutions 


DEPTH(km) 


GENETIC  ALGORITHM  SEARCH 
TEXAS950414  HRV 

S-VEL 

2.0  3.0  4.0  5.0 


GOODNESS  >  0.7  00  BEST=0.919 

Fig.  4.14.  The  GA  searched  models  for  HRV. 


at  a  given  depth  qualitatively  indicates  the  sensitivity  of  the  velocity  to 
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LON=-103.327  LAT=30.261  DEPTH=  23.00 
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Fig.  4.15.  The  waveform  fit  of  the  best  GA  searched  model  for  FRB. 

the  data.  For  example,  for  the  ANMO  transverse  component,  the 
crustal  velocities  are  better  defined  than  upper  mantle  velocities, 
which  is  as  expected  for  a  signal  recorded  only  600  km  away  from  a 
crustal  event. 
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Fig.  4.16.  The  GA  searched  models  for  FRB. 

The  experience  we  learned  here  can  he  used  to  build  a  more 
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Fig.  4.17.  The  waveform  fit  of  the  best  GA  searched  model  for  LMN. 

intelligent  algorithm  for  modeling  surface-wave  waveform  in  the 
future.  This  future  algorithm  should  use  both  amplitude  and  time- 
shift  information  as  its  criteria  of  goodness-of-fit.  Also,  the  better  way 
of  parameterization  will  further  improve  GA  search  technique’s  appli¬ 
cability. 
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Fig.  4.18.  The  GA  searched  model  for  LMN. 
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Fig.  4.19.  The  waveform  fit  for  INK.  INK  is  located  in  the  most  northwest  of  Canada. 
From  the  Rayleigh  wave  trace,  we  can  see  the  wavetrain  due  to  the  multipathing 
effect  come  in  after  the  1600  seconds.  The  current  used  criteria  of  goodness-of-fit  can 
not  properly  represent  the  crust  wave  and  cause  the  wrong  long  lasting  S3mthetic 


wave. 
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Fig.  4.20.  The  GA  searched  models  for  INK.  The  unreasonable  low  crust  velocity 
value  is  causing  by  the  criteria  of  goodness-of-fit  which  can  not  properly  represent 
the  crust  wave. 


CHAPTER  5 

COMPARISON  OF  TECHNIQUES 


In  this  chapter,  four  different  techniques  will  be  used  to  modeling 
the  surface-waveform  which  recorded  at  CCM  (Cathedral  Cave,  Mis¬ 
souri).  These  techniques  are  Genetic  Algorithm,  traditional  surface 
inversion  of  dispersion  data,  linear  waveform  inversion  (Gomberg  and 
Masters,  1988),  and  GSDF  inversion.  The  comparison  of  different 
results  will  illustrate  the  strength  and  weakness  of  each  method.  A 
discussion  will  address  on  what  new  knowledge  can  we  learn  through 
the  waveform  modeling  and  what  is  the  good  way  to  judge  the  reliabil¬ 
ity  of  inversion  results. 

The  reason  to  use  CCM  data  in  this  test  is  that  we  found  CCM 
waveform  is  very  difficult  to  model.  It  can  not  be  fit  by  any  simple 
model.  Through  all  these  waveform  modeling  tests,  we  believe  that 
some  new  knowledge  can  be  learned.  All  these  experiences  can  be 
used  to  improve  the  algorithms  of  waveform  modeling. 

5.1  Genetic  Algorithms 

In  applying  GA  search  method  to  CCM,  two  steps  search  were  per¬ 
formed.  As  mentioned  in  chapter  4,  there  is  a  cycle-skipping  problem 
in  appl3dng  many  layers  GA  search  to  long  distance  stations.  Base  on 
this  reason,  the  first  step  search  were  performed  using  4-layers  GA 
search.  The  best  result  from  this  first  step  search  is  used  to  recon¬ 
struct  the  searching  bounds  for  the  second  step  GA  search. 

Figure  5.1  shown  the  first  step  4-layers  GA  searched  results  and 
searching  bounds.  The  best  model  from  first  step  search  can  match 
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Fig.  5.1.  First  step  4-layers  GA  searched  result. 


the  low  frequency  waveform  (Figure  5.2).  According  the  first  step 
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Fig.  5.2.  The  waveform  fit  for  the  best  model  in  first  step  GA  search. 


search  result,  the  second  step  GA  searching  bounds  were  designed. 
The  second  step  GA  search  result  are  shown  as  Figure  5.3.  The  wave¬ 
form  predicted  from  the  best  model  can  fit  the  fundamental  mode  as 
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Fig.  5.3.  Second  GA  searched  results. 


high  as  0.1  Hz,  but  the  higher  mode  is  out  of  phase.. 
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Fig.  5.4.  The  waveform  fit  of  the  best  model  in  second  step  GA  search. 


The  GA  search  method  provides  a  very  good  starting  model  which 
has  the  fundamental  mode  matched  in  phase.  The  higher  mode  may 
be  out  of  phase,  but  the  model  is  a  good  starting  model  for  final  tuning. 
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The  most  important  advantage  for  GA  search  method  is  that  requires 
no  human  judgement  during  the  search.  The  consumed  time  is  the 
least  among  all  methods  tested  in  this  study. 


RAYLEIGH 


PERIOD 


Fig.  5.7.  The  comparison  of  observed  (circle)  and  predicted  (line)  phase  velocity  of 
inversion. 


5.2  Inversion  of  Dispersion  Data 

The  traditional  siuface  wave  inversion  method  to  get  seismic 
structure  is  using  the  dispersion  data  (e.g.  Russell,  1987).  This 
method  requires  first  to  evaluate  surface  wave  dispersion  data.  The 
Multiple  Filter  Technique  (Dziewonski  et.  al.,  1969;  Herrmann,  1973) 
is  often  used  in  extracting  surface  wave  group  velocity  dispersion 
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Fig.  5.5.  The  extracted  surface  wave  dispersion  data  using  multiple  filtering  tech¬ 
niques.  The  top  plot  shows  the  spectrum  amplitudes  at  different  periods.  The  bottom 
plot  shows  the  group  velocity  dispersion  curve. 
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Fig.  5.6.  This  figure  shows  the  extracted  traces  using  phase  match  filter  technique. 
The  top  trace  is  the  vertical  component  displacement  type  seismogram  of  COM.  The 
second  trace  is  extracted  fundamental  mode.  The  third  trace  is  the  extracted  first 
higher  mode.  The  bottom  trace  is  the  residual. 

curves.  However,  inversion  using  only  group  velocity  dispersion  data 
usually  suffer  nonuniqueness  and  unstability.  The  better  way  is  using 
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Fig.  5.8.  The  inversion  model. 

both  group  and  phase  velocity  dispersion  in  inversion.  This  require 
another  procedure  to  extracting  surface  wave  phase  velocity  dispersion 
data.  Herrin  and  Goforth  (1977)  introduce  the  Phase-Matched  Filter 
method  to  estimate  the  phase  velocity.  After  these  two  procedures,  we 
can  perform  the  traditional  surface  wave  dispersion  curve  inversion. 

In  figure  5.5,  we  show  the  multiple  filtering  result.  The  bottom 
figure  shows  the  group  velocity  dispersion  curve.  In  COM,  the  vertical 
component  of  displacement  type  record  shows  a  very  clear 
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Fig.  5.9a.  The  waveform  fit  of  observed  and  predicted  seismograms.  The  model  is 
shown  as  figure  5.8.  (a)  Displayed  at  frequency  range  0.01-0.1  Hz. 


fundamental  mode  and  first  higher  mode  Rayleigh  wave.  For  period 
longer  than  40  seconds,  the  fundamental  mode  curve  is  not  smooth  as 


105 


TEXAS950414 


CCM  Disp 


Z 

o  b  s  .  Z 
5.169E-03 

s  y  n  .  Z 
8.660E-03 


R 

o  b  s  .  R 
5.577E-03 

s  y  n  .  R 
5.297E-03 


T 

o  b  s  .  T 
9.040E-03 


FILTER  Flo=0.010  Fhi  =  0.500  (Hz)  Norder=  4 
LON=-103.327  LAT=30.261  DEPTH=  23.00 

AZ=48.918  BAZ=235.704  DIST=  1408.160 


Fig.  5.9b.  (Cont’d)  (b)  The  waveform  fit  for  the  frequency  range  0.01-0.5  Hz. 


shorter  periods.  This  character,  which  is  observable  for  stations  in 
similar  geological  region  such  as  CBKS  (shown  on  figure  5.24),  may 
caused  by  some  structure  in  deep.  The  first  higher  mode  is  clearly 


106 


identified  between  6  and  18  seconds.  But  for  periods  shorter  than  5 
seconds,  the  possible  second  higher  mode  is  coming  in  and  form  a 
strong  Lg  phase  at  3.4  km/sec.  We  would  like  to  see  whether  all  these 
features  can  be  modeled  through  waveform  modeling.  Figure  5.6 
shows  the  extracted  fundamental  and  first  higher  mode  using  phase- 
matched  filter  technique.  The  top  trace  is  the  observed  data  and  the 
bottom  trace  is  the  residual  seismogram  after  subtracting  the 
extracted  fundamental  and  first  higher  mode  Rayleigh  wave.  In  the 
residual  seismogram,  we  can  see  a  possible  second  higher  mode  signal 
in  the  time  window  320  to  430  seconds.  But  there  still  remain  one 
question  about  what  kind  signal  is  that  which  arrives  between  440-540 
seconds  with  strong  energy. 

Using  group  and  phase  velocities  in  inversion,  the  results  are 
shown  on  figure  5.7  and  5.8.  Figure  5.7  shows  the  used  data  and  pre¬ 
dicted  curve  by  inverted  model.  Figure  5.8  shows  the  model  of  inver¬ 
sion.  The  starting  model  is  the  best  model  obtained  in  first-step  GA 
search  shown  on  figure  5.1.  Figure  5.9  shows  the  waveform  fit  using 
the  inverted  model  at  two  different  frequency  bands  0.01-0.1  and 
0.01-0.5  Hz.  There  is  a  priori  information  used  in  generated  synthet¬ 
ics  which  is  not  available  from  this  inversion.  That  is  the  frequency 
dependent  Q  which  will  be  described  later  in  this  chapter.  From  the 
dispersion  curve  fit,  it  is  not  surprising  to  see  that  synthetic  funda¬ 
mental  mode  is  slightly  slow  than  observed  for  the  crustal  wave  and 
the  mantle  wave  is  not  matched  (figure  5.9a).  Also,  the  first  higher 
mode  (figure  5.9b)  should  be  fitted  much  better  than  fundamental 
mode. 

This  exercise  shows  that  traditional  surface  wave  inversion  can 
provide  a  fairly  good  model  from  the  point  of  view  of  waveform  model¬ 
ing.  However  this  method  also  has  its  limitations  such  as  the  inver¬ 
sion  has  been  dominated  by  short  period  data  which  has  more  data 
point  than  long  period  data  and  has  no  ability  to  change  model  to  fit 
the  data  around  30  seconds  in  this  case.  However,  this  technique  is 
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very  fast  comparing  to  waveform  inversion  techniques. 
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Fig.  5.10.  The  waveform  comparison  of  observed,  sjmthetic,  and  the  residual  seismo¬ 
grams  at  two  different  frequency  bands.  The  top  panel  which  corresponding  to  appar¬ 
ent  velocity  3. 3-2.5  km/sec  is  plotted  at  0.02-0.05  Hz.  The  bottom  panel  for  4.5-2.5 
km/sec  at  0.02-0.1  Hz. 


5.3  Linear  Waveform  Inversion 

Gomberg  and  Masters  (1988)  introduce  the  linear  waveform  inver¬ 
sion  technique.  This  technique  formulates  the  residual  seismogram  of 
observed  and  synthetic  seismograms  in  terms  of  model  parameter  per¬ 
turbations.  In  this  section,  we  will  show  the  inversion  result  obtained 
by  appl5dng  the  method  to  CCM  data.  The  programs  used  here  is 
coded  by  professor  Robert  Herrmann,  Saint  Louis  University. 

In  the  original  paper  (Gomberg  and  Masters,  1988),  the  inversions 
were  performed  directly  using  whole  frequency  range  interested. 
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However,  it  is  found  by  professor  Herrmann,  that  a  better  way  to  using 
this  tool  in  inversion  is  performing  the  inversions  in  different  fre¬ 
quency  bands.  This  improved  procedure  can  utilize  different  informa¬ 
tion  for  low  frequency  and  high  frequency  to  constrain  the  model  in  a 
more  stable  and  robust  way.  The  similar  idea  was  reported  by  Bimks 
et  al.  (1995)  for  studying  the  reflection  data. 

Using  the  best  model  obtained  from  first-step  GA  search,  and  per¬ 
form  the  linear  waveform  inversion.  In  figure  5.10,  it  shows  the  wave¬ 
form  fit  at  two  frequency  ranges.  We  can  see  that  the  phase  and 
amplitude  of  fundamental  mode  is  matching  the  observed.  The  first 
higher  mode  does  not  match  the  observed  data  because  the  a  priori 
information  of  frequency  dependent  Q  is  not  used  in  this  test.  The 
inverted  model  is  not  shown  because  it  is  used  to  do  fine  adjustment  by 
trial  and  error  and  a  better  waveform  fit  is  obtained  which  will  be  pre¬ 
sented  later. 

Although  all  reported  applications  of  linear  waveform  inversion 
were  for  the  short  distance  records,  we  believe  it  is  possible  to  apply 
this  technique  to  longer  distance  records.  The  strength  for  this 
method  is  the  inversion  will  force  on  the  large  residual  between 
observed  and  predicted  data.  However,  the  weakness  of  this  method  is 
associated  with  its  strength,  i.e.  it  needs  a  sophisticated  controlling 
method  to  avoid  the  effects  that  caused  by  improper  knowledge  of 
source. 

5.4  GSDF  Inversion 

In  this  section,  we  are  appl3dng  the  GSDF  inversion  method  to 
CCM  data.  As  the  model  obtained  from  linear  waveform  inversion,  the 
inverted  model  is  not  shown  here,  and  an  improved  waveform  fit,  by 
trial  and  error  adjustment,  based  on  this  model  will  be  shown  later  in 
this  chapter.  Figure  5.11  shows  the  waveform  fit  of  the  GSDF  inver¬ 
sion  result. 
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Fig.  5.11a.  The  waveform  fit  for  the  model  obtained  from  GSDF  inversion.  Three  dif¬ 
ferent  frequency  ranges  are  shown:  (a)  for  0.01-0.05  Hz. 
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Fig.  5.11b.  (Cont’d)  (b)  for  0.01-0.1  Hz. 


5.5  How  To  Judge  the  Waveform  Fit 
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Fig.  5.11c.  (Cont’d)  (c)  for  0.01-0.5  Hz. 


What  is  the  criteria  to  assess  the  goodness  of  waveform  fit?  From 
the  viewpoint  of  waveform  modeling,  the  synthetics  should  be  able  to 
match  the  observed  data  in  both  phase  and  amplitude  for  different 
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Fig.  5.12.  The  partial  derivatives  of  first  higher  mode  at  10  seconds  which  computed 
for  GSDF  inversion. 


frequency  bands.  In  this  case,  first  let’s  see  what  are  observable  from 
recorded  seismograms  of  station  CCM.  In  CCM,  we  can  see  fundamen¬ 
tal  mode  and  first  higher  mode  with  large  amplitude,  so  the  first 
requirement  should  be  matching  both  fundamental  and  higher  modes 
signal.  If  we  can  find  such  a  model  which  can  predict  observed  surface 
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Fig.  5.13a.  Waveform  fit  for  the  model  which  increasing  the  upper  mantle  Q.  Three 
different  frequency  ranges  are  displayed:  (a)  0.01-0.05  Hz. 


waveform,  we  would  like  to  know  whether  this  model  is  good  enough  to 
match  the  body  wave,  P^i  and  S  in  this  case.  After  that,  we  would  like 
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Fig.  5.13b.  (Cont’d)  (b)  0.01-0.1  Hz. 


to  know  if  a  model  can  generally  meet  the  first  two  requirement,  is 
there  a  way  to  tell  us  how  confident  this  model  can  be  or  any  improve¬ 
ment  can  be  expected  in  the  future  modeling. 


115 


TEXAS950414 


CCM  Disp 


Z 

o  b  s  .  Z 
5.711E-03 

S  V  11  .  Z 
9.''630E-03 


R 

o  b  s  .  R 
5.551E-03 

s  y  n  .  R 
6.303E-03 


T 

o  b  s  .  T 
8.588E-03 


FILTER  Flo=0.010  Fhi=0.500  (Hz)  Norder=  4 

LON=-103.327  LAT=30.261  DEPTH=  23.00 


AZ=48.918  BAZ=235.704  DIST=  1408.160 


Fig.  5.13c.  (Cont’d)  (c)  0.01-0.5  Hz. 


From  above,  we  can  say  that  all  four  algorithms  can  provide  fairly 
good  models  that  can  match  fundamental  mode  up  to  0.1  Hz.  How¬ 
ever,  those  model  can  not  fit  the  higher  mode,  therefore  the  question  is 
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Fig.  5.14.  This  figure  shows  a  better  model  which  based  on  GSDF  inversion  result, 
and  using  trial  and  error  to  do  the  fine  timing. 

how  to  match  the  higher  mode.  Especially  looking  at  the  GSDF  inver¬ 
sion  result  (figure  5.11b),  we  can  see  predicted  higher  mode  amplitude 
is  about  factor  two  smaller  than  observed  data. 

There  are  two  possible  ways  to  increase  the  higher  mode  ampli¬ 
tude.  Let’s  check  all  the  partial  derivatives  created  for  the  GSDF 
inversion.  As  shown  on  figure  5.12,  we  can  find  that  if  we  increase  the 
upper  mantle  Q  value,  it  is  possible  to  increase  the  higher  mode  ampli¬ 
tude.  Of  course,  the  crustal  Q  value  should  decrease  to  compensate 
the  change  and  to  maintain  the  fundamental  mode  amplitude.  The 
result  is  shown  as  figure  5.13.  Indeed,  the  higher  mode  amplitude  is 
increasing,  but  compare  to  figure  5.11  we  can  find  not  only  higher 
mode  but  also  fundamental  mode  amplitude  are  increasing.  Especially 
the  low  frequency  part  of  fundamental  mode  which  corresponding  to 
the  wave  traveling  through  the  upper  mantle.  Another  problem  is 
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Fig.  5.15b.  (Cont’d)  (b)  0.01-0.1  Hz. 


The  second  way  to  increase  the  higher  mode  amplitude  is  to  con¬ 
sider  frequency  dependent  Q  for  high  frequency  signal.  To  preserve 
the  low  frequency  fundamental  mode  amplitude,  the  frequency 
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FILTER  Flo=0.010  Fhi  =  0.500  (Hz)  Norder=  4 

LON=-103.327  LAT=30.261  DEPTH=  23.00 

AZ=48.918  BAZ=235.704  DIST=  1408.160 

Fig.  5.15c.  (Cont’d)  (c)  0.01-0.5  Hz. 

dependent  Q  only  apply  to  signal  whose  frequency  content  is  greater 
than  0.1  Hz.  To  avoid  the  abrupt  change  in  spectrum  shape  when 
appl5dng  frequency  dependent  Q,  a  transition  zone  is  used  for  0.1-0.2 
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Fig.  5.17.  This  figure  shows  a  better  model  which  based  on  linear  waveform  inversion 
result,  and  using  trial  and  error  to  do  the  fine  tuning. 

Hz.  The  following  two  examples  are  based  on  trial  and  error  adjust¬ 
ment  for  two  models  from  GSDF  inversion  and  linear  waveform  inver¬ 
sion.  The  reference  period  of  Q  is  0.1  Hz  and  rj  =  0.5.  During  the  trial 
and  error  adjustment,  the  fundamental  mode  phase  may  shift  and  out 
of  phase.  But  the  purpose  is  to  see  how  to  fit  the  higher  mode  phase 
and  amplitude,  and  what  process  creates  such  strong  Lg  wave. 

The  first  example  is  a  fine  tuning  result  for  the  GSDF  inversion 
result.  The  model  is  shown  as  figure  5.14  and  using  fi'equency  depen¬ 
dent  Q.  The  waveform  fit  is  displayed  as  figure  5.15.  From  the  wave¬ 
form,  we  can  see  that  low  frequency  fundamental  mode  amplitude  does 
not  change  much.  Through  the  trial  and  error  adjusting,  the  higher 
mode  phase  and  amplitude  are  generally  matched.  However,  the  Airy 
phase  amplitude  is  increased  during  such  adjustment.  The  waveform, 
again,  can  not  match  the  signals  in  the  time  window  420-440  seconds; 
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Fig.  5.18a.  The  waveform  fit  of  the  fine-tuning  model  shown  at  figure  5.17.  Three  dif¬ 
ferent  frequency  ranges  are  shown:  (a)  0.01-0.05  Hz. 


also  it  can  not  explain  what  high  frequency  signals  are  those  between 
440  and  470  seconds.  Looking  at  the  single  modes  display  (figure 
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Fig.  5.18c.  (Cont’d)  (c)  0.01-0.5  Hz. 

The  second  example  is  based  on  the  inversion  result  of  linear 
waveform  inversion.  This  adjusted  model  shows  good  agreement  with 
observed  data,  so  we  will  use  this  model  as  an  example  to  illustrate 
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Fig.  5.19.  The  waveform  fit  which  computed  using  locked  mode  approximation  for 
the  fine-tuning  model  shown  at  figure  5.17. 


what  is  the  reliability  of  model.  In  this  model,  not  only  the  Rayleigh 
wave  but  also  the  Love  wave  are  modeled.  Figure  5.17  shows  the 
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Fig.  5.20.  The  observed,  synthetic,  and  seven  single-mode  seismograms  were  shown 
here. 


models  for  the  Rayleigh  and  Love  waves.  There  are  several  characters 
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Fig.  5.21.  The  resolving  kernel  for  the  model  shown  on  figure  5.18. 

worth  to  mention.  First  of  all,  during  the  adjusting  process,  no  low 
velocity  zone  is  allow  to  exist  in  the  crust.  Second  of  all,  there  is  no 
need  to  introduce  a  funny  low  velocity  zone  in  the  uppermost  mantle 
(40-220  km)  which  is  usually  seen  in  other  reported  surface  wave 
inversion  results.  Third  of  all,  the  model  shows  a  gradient  zone 
between  40  and  60  km  instead  of  a  sharp  crust-mantle  boundary. 
Final,  the  lower  crust  (20-40  km)  is  a  steep  velocity  zone.  Figure  5.18 
shows  the  waveform  fit  at  three  different  fi^equency  ranges.  From 
these  plots,  we  can  find  that  the  higher  mode  waveform  fit  the 
observed  data  very  well.  Especially  the  signals  in  the  400-440  time 
window  which  is  not  well  modeled  by  the  other  models.  From  figure 
5.20,  the  single  modes  display  shows  a  very  good  first  higher  model 
envelop  with  a  nice  tail  after  410  seconds.  This  is  consistent  with  the 
observed  data.  It  also  show  that  the  amplitude  variation  is  possibly 
caused  by  construction  and  destruction  of  higher  modes.  Figure  5.19 
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Fig.  5.22.  The  measurements  of  four  GSDF  functionals  for  the  model  shown  on  figure 
5.18. 


shows  the  synthetic  seismograms  which  generated  by  locked  mode 
approximation  can  match  the  P^i  phase.  One  of  the  importances  of 
waveform  modeling  shown  here  is  that  it  demonstrates  a  single  1-D 
model  can  match  both  body  wave  and  surface  wave  waveforms  in  a 
regional  stable  area. 

How  confident  can  this  model  tell?  In  this  study,  the  inversions 
were  performed  using  damped  least  square  method  (Marquardt,  1963). 
The  calculation  of  resolving  kernels  are  straightforward.  Figure  5.21 
shows  the  resolving  kernels  and  standard  deviation  of  model  parame¬ 
ters.  However,  the  standard  deviation  bar  is  depend  on  the  damping 
value  used,  it  really  can  not  show  any  information.  On  the  other  hand. 
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Fig.  5.24.  The  waveform  fit  for  station  CBKS  at  frequency  band  0.01-0.5  Hz. 

the  measxirements  of  four  GSDFs  may  provide  extra  information  to 
indicate  the  goodness  of  waveform  fit.  Figure  5.22  show  the  four 
GSDF  measurements  of  vertical  component  Rayleigh  wave.  From  the 
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Fig.  5.25.  The  fine-tuning  model  for  station  CBKS  which  based  on  the  GSDF  inver¬ 
sion  result. 

group  velocity  differential  time  dTG,  it  indicate  that  around  37  seconds 
which  represents  the  predicted  Airy  phase  envelop  is  15  seconds  slow 
than  observed.  From  dTP,  it  is  clear  to  see  that  this  model  may  doing 
well  in  high  frequency  but  the  low  frequency  part  is  systematically  5 
seconds  faster  than  data.  From  dTQ,  we  can  see  that  for  the  period 
around  35  seconds,  the  synthetic  seismogram’s  amplitude  is  too  large. 
The  GSDF  measurements  are  better  indicator  of  waveform  fit  than 
resolving  kernels. 

So  far,  there  are  two  new  hypotheses  proposed.  The  first  hypothe¬ 
sis  is  that  a  gradient  velocity  transition  zone  exist  at  the  uppermost 
mantle  instead  of  a  sharp  crust-mantle  boundary.  The  second  hypoth¬ 
esis  is  that  whether  the  frequency  dependent  Q  is  real  or  not.  We  will 
use  another  two  stations  in  the  same  region  to  test  these  two  h3rpothe- 


ses. 
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Fig.  5.27.  The  waveform  fit  at  0.01-0.5  Hz  of  station  WMOK  for  the  model  shown  on 
figure  5.28. 


The  first  station  is  CBKS.  From  its  group  velocity  dispersion 
curve  (figure  5.23),  it  shows  almost  no  interference  from  noise.  The 
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Fig.  5.28.  The  fine-timing  model  for  WMOK  which  based  on  GSDF  inversion  result. 


GSDF  inversion  result  shown  as  figure  5.25  and  its  associated  wave¬ 
form  fit  is  shown  as  figure  5.24.  The  second  station  used  to  test  those 
hypothesis  is  WMOK.  The  hypocentral  distance  is  very  short,  only 
about  650  km.  The  dispersion  curve  is  not  as  clean  as  CBKS,  only  the 
8-30  seconds  fundamental  mode  can  be  identified.  The  waveform  fit  is 
shown  as  figure  5.27  and  the  model  is  shown  as  figure  5.28.  Although 
the  waveform  fit  for  both  CBKS  and  WMOK  are  not  perfect  for  higher 
mode,  but  they  have  correct  amplitudes.  From  the  models,  they  do  not 
require  frequency  dependent  Q  to  fit  the  amplitudes  of  both  fundamen¬ 
tal  and  higher  modes.  But  comparing  to  CCM,  these  two  models  all 
have  a  high  crustal  Q  values.  However,  there  are  two  possibilities  for 
this: 

(1).  The  frequency  dependent  Q  behavior  is  not  observable  for  short 
distance  stations. 
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(2).  The  observed  frequency  dependent  Q  at  CCM  is  an  artifact  due  to 
improper  source  radiation  pattern  which  create  large  amplitude 
synthetics  and  force  the  crustal  Q  decreasing  to  compensate  this 
effect  and  frequency  dependent  Q  behavior  is  then  required  to 
match  the  higher  mode  amplitude. 

Although,  there  is  no  conclusive  answer  for  this  question,  but  this 
question  which  related  to  the  fundamental  behavior  of  earth  structure 
is  noticed  through  modeling  waveforms. 

By  comparing  models  of  CCM,  CBKS,  and  WMOK,  it  is  found  that 
all  three  models  have  a  velocity  transition  zone  at  the  uppermost  man¬ 
tle.  Durrheim  and  Mooney  (1994)  use  seismic  and  geochemical  data  to 
constrain  the  evolution  of  Precambrian  lithosphere  evolution.  Their 
conclusion  is  that  underplating  is  only  happen  at  the  Proterozonic 
instead  Archean  crust-mantle  boundary.  The  seismic  signature  of 
underplating  is  a  basal  layer  with  high  P-velocity  (7.0-7. 6  km/sec).  In 
oxir  study,  the  source  is  located  between  the  Proterozonic  platform  and 
Rocky  mountains,  and  the  station  CCM,  CBKS,  and  WMOK  lay  inside 
the  Proterozonic  platform,  therefore  the  propagation  paths  are  sam¬ 
pling  the  relatively  uniform  platform.  With  fixed  Poisson’s  ratio  in 
inversion,  the  transitional  velocity  zone  between  40  and  60  km  shown 
in  each  model  should  be  the  evidence  of  underplating  process.  It  would 
be  interesting  to  study  the  Archean  region  using  waveform  modeling 
techniques  to  see  whether  a  basal  layer  exists  or  not. 

There  is  one  question  remain  unanswered.  That  is  what  are  those 
signals  arrives  between  430  and  470  seconds  in  vertical  component 
CCM  data.  We  have  no  answer  for  this.  But  from  figure  5.9(b),  the 
dispersion  curve  inversion  result,  we  see  a  clue. 

From  results  of  reflection  and  refraction  studies  of  continental 
lithosphere,  it  has  been  observed  that  there  exists  a  reflective  lower 
crust  (e.g.  Mooney  and  Brocher,  1987).  Such  seismic  reflectors  can  be 
related  to  geological  interpretations  such  as  thrust  stack,  duplex,  and 
so  on  (Hatcher,  1986).  Braile  and  Chiang  (1986)  conduct  numerical 
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experiments  to  study  what  kind  structure  can  produce  the  reflective 
Moho.  Their  tests  show  that  a  Moho  transition  zone  with  many  fine 
layers  which  have  small  velocity  variations  (ZIG-ZAG  pattern)  can  be 
the  possible  model. 

During  the  trial  and  error  modeling  process,  we  manually  keep 
the  velocity  model  as  smooth  as  possible  and  remove  the  low  velocity 
zones.  The  layer  thickness  of  the  model  is  5  km  for  top  60  km.  For 
such  layer  thickness  and  smooth  velocity  curve,  it  will  only  generate  a 
smooth  fundamental  mode.  Therefore,  a  fine-layer  velocity  model  with 
small  perturbations  is  a  possible  origin  for  those  signals. 

5.6  Discussion 

In  this  chapter,  we  test  four  different  algorithms  to  model  the  sur¬ 
face  waveform.  The  results  turn  out  each  algorithm  can  provide  differ¬ 
ent  waveform  fit  result  which  emphasize  different  part  of  waveform 
due  to  their  different  criteria  and  formulations.  For  example,  the  GA 
search  algorithm  can  find  models  that  match  the  most  energetic  phase, 
fundamental  mode  in  this  case,  but  it  lacks  the  ability  to  match  the 
other  less  energetic  phases  such  as  higher  modes.  The  reason  is  sim¬ 
ple  because  of  the  criteria  used.  In  a  similar  fashion,  the  linear  wave¬ 
form  inversion  is  also  controlled  by  the  large  residuals  which  might  be 
artifacts  caused  by  improper  source.  Using  different  windows  and  fil¬ 
ters  to  force  the  inversion  to  focus  on  different  part  of  waveform  is  a 
possible  solution,  but  the  controlling  mechanism  will  strongly  depend 
on  user’s  expertise.  The  traditional  dispersion  curve  inversion  also 
has  its  own  problems  such  as  the  accuracy  in  phase  velocity  measure¬ 
ments  and  the  data  density  used  in  inversion.  Like  the  traditional  dis¬ 
persion  curve  inversion,  GSDF  inversion  also  suffer  fi'om  data  density 
problem  because  both  of  them  utilize  the  different  remeasurements 
from  waveform. 
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All  four  algorithms  can  provide  reasonable  model  with  waveform 
fit  as  good  as  up  to  0.1  Hz.  However,  more  structural  information  are 
revealed  in  higher  frequency  contents,  the  design  of  a  model  refining 
algorithm  will  be  the  task  for  seismologists.  From  this  study,  we 
believe  that  a  multi-isolation-filters  GSDF  inversion  algorithm  can  be 
a  possible  solution. 

What  is  a  good  criteria  for  waveform  modeling?  So  far,  there  is  no 
answer  for  this  question.  But  we  know  what  is  not  a  good  criteria.  To 
measure  a  simple  pulse,  several  widely  used  criteria  like  cross¬ 
correlation  and  norm-2  can  serve  the  purpose.  However  to  judge  the 
similarity  of  a  complicated  waveform,  all  these  criteria  will  not  give  a 
satisfactory  answer.  The  cross-correlation  can  only  tell  the  similarity 
of  the  most  energetic  phase.  The  2-norm  will  behave  as  a  oscillatory 
wavelet,  like  cross-correlation,  when  the  cycle-skipping  problem 
involved.  Therefore,  we  believe  that  the  success  of  waveform  inversion 
algorithms  will  depend  on  how  the  criteria  have  been  designed. 

When  a  good  model  obtained,  without  tectonic  interpretation  it  is 
just  a  model.  We  would  like  outline  the  possible  contributions  from 
interpreting  waveform  modeling  results.  First  we  would  ask  what  is 
the  physical  process  happen  in  continental  lithosphere?  So  far,  the  evi¬ 
dences  from  all  branchs  of  geoscience  are  too  few  to  give  us  a  conclu¬ 
sive  picture  about  what  happen  down  there.  So  we  will  ask  what  kind 
seismic  evidences  can  we  provide  to  this  topic? 

Jordan  (1975)  proposed  the  continental  tectosphere  hypothesis. 
The  continental  tectosphere  will  have  a  cool  thick  root  (up  to  400  km). 
Using  nmnerical  modeling,  the  influence  of  such  thick  tectosphere  on 
the  geodynamic  process  of  mantle  and  on  the  plate  motion  has  been 
modeled  (e.g.  Stoddard  and  Abbott,  1996).  Therefore,  the  first  ques¬ 
tion  is  how  to  determine  the  thickness  of  continental  lithosphere? 

As  described  previously,  the  evolution  model  of  Precambrian  litho¬ 
sphere  proposed  by  Durrheim  and  Mooney  (1994)  suggest  that  the 
Archean  lithosphere  is  thicker  than  Proterozonic  lithosphere.  Some 
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mantle  process  are  different  such  as  underplating  only  happen  in  Pro- 
terozonic  crust-mantle  boundary.  Therefore,  the  second  question  will 
be:  Can  we  find  the  fundamental  difference  of  these  two  regions  from 
seismological  studies? 

Sato  et  al.  (1989)  suggest  that  seismic  velocity  and  Q  model  can  be 
used  to  estimate  the  temperature  gradient  and  partial  melt  fraction  in 
the  upper  mantle.  The  temperature  gradient  can  be  inferred  from  the 
waveform  modeling  results.  This  could  be  another  indirect  constraint 
besides  the  heat  flow  measurements  for  continental  hthosphere  evolu¬ 
tion. 

Anisotropy  is  a  common  observed  phenomenon  from  waveform 
modeling  result.  From  observing  the  SKS  and  SKKS  splitting,  Silver 
and  Chan  (1991)  suggest  that  the  anisotropy  is  caused  by  the  subconti¬ 
nental  upper  mantle  deformations  during  the  different  tectonic 
episodes.  Therefore,  anisotropy  can  be  another  indicator  to  the  evolu¬ 
tion  history  of  continental  lithosphere. 


CHAPTER  6 

DISCUSSION  AND  CONCLUSION 


In  this  study,  several  algorithms  have  been  used  to  model  surface 
waveform.  Each  algorithm  has  different  advantages  and  weaknesses. 
One  of  the  most  important  conclusion  is  that,  through  carefully  wave¬ 
form  modeling,  both  body  wave  and  surface  wave  can  be  modeled  using 
1-D  model.  Waveform  modeling  can  also  provide  more  detail  structure 
information  which  is  not  obtainable  by  fitting  gross  data  sets  such  as 
surface  wave  dispersion  data. 

The  difference  of  tested  algorithms  are  mainly  in  the  "data"  prepa¬ 
ration  for  inversion.  The  linear  inversion  scheme  used  is  working  for 
inversion  with  no  a  priori  information  (Lawson  and  Hanson,  1974). 
Some  issues  that  involving  the  different  inversion  schemes  such  as 
incorporating  a  priori  information  (Tarantola  and  Valette,  1982),  and 
constrained  inversion  (Carrion,  1989)  has  not  been  tested  in  this  study. 

Another  issue  that  related  to  inversion  schemes  is  the  weighting 
functions.  For  linear  inversion,  each  parameter  is  acting  as  a  free 
parameter,  usually  the  inversion  will  become  unstable  due  to  the 
parameterization  and  data  density  used  in  inversion.  Russell  (1987) 
introduced  "differential"  weighting  function  to  limit  the  parameters 
movement.  However,  when  dealing  with  teleseismic  records  or  many 
thin  layers  parameterization,  this  weighting  function  will  not  be  able 
to  avoid  producing  some  funny  low-velocity  zones  in  model.  Therefore, 
a  Nolet  style  weighting  function  (Nolet,  1990;  Snoke  et  al. ,  1997)  could 
be  tested  in  the  future  to  see  can  it  stablize  the  inversion  result. 

what  can  be  done  in  future?  In  the  aspect  of  technique  improve¬ 
ment:  First  of  all,  from  this  study,  we  notice  that  the  envelope  of 
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surface  wave  is  a  good  indicator  for  Q  structure.  We  believe  that  to 
formulate  an  algorithm  to  invert  Q  structure  from  different  band- 
passed  surface  wave  envelop  will  be  a  better  tool  to  study  upper  man¬ 
tle  Q  structure.  The  similar  idea  has  already  been  reported  by  Cara  et 
al.  (1987)  and  Nolet  (1990).  Second  of  all,  a  fine  tune-up  algorithm  is 
needed  to  match  the  less  energetic  higher  mode  waveforms  without 
destroy  the  matched  fundamental  mode.  We  believe  it  can  be  imple¬ 
mented  by  constructing  a  multiple  isolation  filters  GSDF  inversion 
scheme  which  will  utilize  information  from  different  wavegroups 
simultaneously.  Another  alternative  will  be  multi  windowed  linear 
waveform  inversion.  Third  of  all,  to  implement  the  anisotropy  inver¬ 
sion  for  surface  waveform.  Finally,  to  invent  a  inversion  using  both 
both  surface  and  body  waves. 

In  the  aspect  of  application:  We  strongly  feel  that  waveform  inver¬ 
sion  techniques  can  be  very  powerful  tools  to  study  regional  events. 
First  of  all,  regional  events  have  their  energy  concentrate  in  a  short 
time  window  which  may  have  less  influence  by  multi-pathing  effects. 
The  GA  search  results  shown  in  chapter  4  can  support  this  argument. 
Second  of  all,  the  propagation  path  is  more  likely  travel  through  the 
same  geological  region.  Although  such  modeling  may  only  provide 
information  of  average  crustal  structure,  but  it  can  apply  to  many 
regions.  It  would  be  interesting  to  apply  waveform  inversion  methods 
in  a  more  small  region. 

In  the  aspect  of  global  search  technique:  The  major  thing  need  to 
be  done  is  to  find  a  better  criteria  of  waveform  fit.  This  criteria  is 
directly  concern  the  success  of  global  search  technique  for  waveform 
modeling.  The  second  work  is  to  find  a  way  to  combine  GA  search  and 
gradient  information  to  speed  the  process.  This  will  increase  the 
application  of  global  search  technique.  The  other  thing  is  to  test  other 
parameterization  for  GA  search. 
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